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ABSTRACT 


A finite element program has been developed to perform quasi-static, elastic-plastic crack 
growth simulations. The model provides a general framework for mixed-mode I/II elas- 
tic-plastic fracture analysis using small strain assumptions and plane stress, plane strain, 
and axisymmetric finite elements. 

Cracks are modeled explicitly in the mesh. As the cracks propagate, automatic remesh- 
ing algorithms delete the mesh local to the crack tip, extend the crack, and build a new 
mesh around the new tip. State variable mapping algorithms transfer stresses and dis- 
placements from the old mesh to the new mesh. The von Mises material model is im- 
plemented in the context of a non-linear Newton solution scheme. The fracture criterion 
is the critical crack tip opening displacement, and crack direction is predicted by the 
maximum tensile stress criterion at the crack tip. The implementation can accommodate 
multiple curving and interacting cracks. An additional fracture algorithm based on nodal 
release can be used to simulate fracture along a horizontal plane of symmetry. A core of 
plane strain elements can be used with the nodal release algorithm to simulate the triaxial 
state of stress near the crack tip. 

Verification and validation studies compare analysis results with experimental data and 
published three-dimensional analysis results. Fracture predictions using nodal release for 
compact tension, middle-crack tension, and multi-site damage test specimens produced 
accurate results for residual strength and link-up loads. Curving crack predictions using 
remeshing/mapping were compared with experimental data for an Arcan mixed-mode 
specimen. Loading angles from 0 degrees to 90 degrees were analyzed. The maximum 
tensile stress criterion was able to predict the crack direction and path for all loading an- 
gles in which the material failed in tension. Residual strength was also accurately pre- 
dicted for these cases. 
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Chapter 1 


Introduction 


Cracks don’t always grow straight ahead; sometimes they turn. Researchers in the frac- 
ture field have known this since they started observing and working to understand the 
fracture behavior of engineering materials. But the means to predict the behavior of 
multiple turning cracks under complex loading situations in the presence of significant 
plasticity have not been previously available. 

Ever since the implementation of the finite element method (FEM) on computers, engi- 
neers have applied the method to increasingly complex problems. Indeed, more than 25 
years ago FEM was first used to simulate elastic-plastic fracture (Anderson, 1973). 
However, applications of FEM to elastic-plastic mixed-mode t fracture problems have 
been limited. There are three fundamental issues that need to be resolved: 1) a fracture 
criterion applicable under mixed-mode loading conditions; 2) a direction criterion to pre- 
dict the direction of growth; 3) a means of representing a growing crack without knowing 
a priori where the crack will go. In the case of linear elastic fracture mechanics (LEFM), 
solutions to each of these issues have already been demonstrated. For elastic-plastic 
fracture mechanics (EPFM), until recently, solutions to the first two issues had not been 
experimentally verified and solutions to the third issue have not been addressed ade- 
quately. This work addresses the third issue: representation that allows curvilinear crack 
growth without the crack path known a priori. 

1.1 Background 

Implicit in any discussion of fracture is the assumed relation of fracture to failure. While 
in some cases fracture is planned (such as when the pull-top of a beverage can is acti- 
vated), many times fracture is not desirable, and contingency plans are in place to prevent 
fracture. However, fracture and its relationship to failure will always exist, even though 
the details of the relationship will change over time as technology evolves. 


t 


Unless explicitly noted, the term mixed-mode refers to Mode-I/II. 
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One guiding force for fracture research during the past decade has been the NASA Air- 
craft Structural Integrity Program (Harris, 1990), but in fracture terminology, the initia- 
tion event was a failure. In the spring of 1988 Aloha Airlines Flight 243 suffered major 
structural failure during flight over Hawaii. The pilot was able to land the plane safely 
with minimal loss of life or injury to the passengers and crew, but the incident caused 
great concern over the safety of the aging commercial aircraft fleet. This concern 
prompted the U. S. Government to fund research at both NASA and the FAA to investi- 
gate the issues involved with understanding the problem of airframe structural integrity, 
and to formulate a plan of action to transfer the knowledge gained to industry. The re- 
sulting NASA program, under which this research is funded, is the NASA Airframe 
Structural Integrity Program at NASA Langley Research Center (Harris, 1990). 

1.2 NASA Analysis Methodology 

The NASA Airframe Structural Integrity Program (NASIP) is designed to gain under- 
standing of the structural integrity problem in many areas, which are broken down into 
two major elements: (a) Quantitative Nondestructive Inspection Technology and (b) 
Structural Integrity Analysis Methodology (Harris, 1996). 

The overall goal of the Quantitative Nondestructive Inspection Technology (NDI) ele- 
ment of the program is to develop nondestructive inspection technology to detect dis- 
bonding, corrosion, and fatigue cracks in fuselage structure (Harris, 1996). The current 
research is not part of the NDI element of the NASIP program, so NDI will not be dis- 
cussed further. 

The overall goal of the Structural Integrity Analysis Methodology is defined from two 
different perspectives: that of materials characterization, whose goal is to develop fatigue 
crack growth and fracture criteria for thin-sheet riveted metallic structures, and that of 
structural analysis, whose goal is to develop and verify the structural analysis methodol- 
ogy for stiffened shells with damage (Harris, 1996). The current research is most closely 
tied to the structural analysis methodology, since our task is to develop software that will 
become part of the structural analysis methodology. However, the work is also closely 
tied to the materials characterization goals of the NASIP program, since the analysis 
methodology implemented must accurately represent the current understanding of how 
the materials actually behave, not only in the laboratory at the coupon and component 
levels, but also at the structural system level. 
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1.3 FRANC2D/L - An Engineering Tool 

There are many deliverables that will result from the NASIP research project. One class 
of deliverables is to have a range of engineering tools in the form of computer codes 
available for engineers in research institutions and industry. Every computer code is not 
intended to solve all problems; each code is intended to focus on a specific range of 
problems. The current research effort at KSU is centered around the FRacture ANalysis 
Code-2D with Layers (FRANC2D/L) computer code. This software was originally de- 
veloped at Cornell University (Wawrzynek and Ingraffea, 1987). 

The Cornell version of FRANC solves two-dimensional linear elastic fracture mechanics 
problems using the finite element approximation. The FRANC code was created in re- 
sponse to a need to model the true evolutionary behavior of crack propagation using finite 
elements, while addressing the limitations of the then current approach to modeling such 
problems (Saouma and Ingraffea, 1981; Swenson, 1985; Gerstle, 1986). The author of 
the code took a novel approach to implementing the finite element database. Rather than 
save a traditional finite element database (simple nodal coordinates and element connec- 
tivity), an advanced topological database is used to save the relationships between all ad- 
jacent nodes, elements and element edges. The result is that during the crack propagation 
phase of an analysis, where new surface area is created in the finite element mesh, the 
topology of the mesh is well defined. This is true even for intermediate operations such 
as local mesh editing, where the mesh is deleted in the region around the crack tip, the tip 
is extended, and a new mesh is generated for the updated configuration. The result is that 
during crack propagation the geometry of the problem is preserved automatically, and 
local remeshing operations on the region around the crack face can operate directly on the 
data base (rather than searching for necessary data) since all of the topology of the prob- 
lem is preserved. 

The current FRANC2D/L version of the software is a direct extension of the Cornell ver- 
sion of the software and is the result of two recent research projects at KSU. Sudhir 
Gondhalekar extended the original two-dimensional code to model layered structures 
such as lap joints and adhesively bonded joints (Gondhalekar, 1992). The model was 
limited to in-plane loading and deformations only. Srinivas Krishnan subsequently ex- 
tended the code to include out-of-plane loading and deformations through linear, small 
strain bending (Krishnan, 1994). The current version of the code retains all of the neces- 
sary pre- and post-processing functionality from the Cornell version, as well as all of the 
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fundamental fracture mechanics functionality, in addition to the new functionality added 
by Gondhalekar and Krishnan. 


1.4 Crack Growth Methodology 

Figure 1.1 shows the conceptual crack growth methodology. This methodology is an 
augmented version of that presented by Potyondy (1993). In the current version pre- and 
post-procession are included for clarity, and the crack growth model will be elaborated 
upon in greater detail. 

The general methodology starts with the pre-processing stage, where the geometry, mesh, 
material properties, and boundary conditions are specified. The result of preprocessing is 
the initial representation, R h In the case of the FRANC2D/L implementation, which uses 
a topological database, this representation for the analysis contains all of the necessary 
information to run an analysis, but organizes the data for efficient pre- and post- 
processing. 

The Analysis Model Generation phase of the methodology extracts data from the data- 
base and converts the data into a form that is efficient for the analysis phase. This phase 
includes such steps as building a traditional node and connectivity list from the database, 
optimizing the equation numbering to minimize the stiffness matrix profile, and initializ- 
ing all of the data structures necessary for the Solution Procedure. 

The Solution Procedure is responsible for applying boundary conditions and satisfying 
the equilibrium equations. The solution procedure is nonlinear incremental-iterative. As 
the solution progresses, the Solution Procedure is also responsible for monitoring fracture 
parameters at appropriate times during the analysis, and responsible for initiating appro- 
priate actions when fracture criteria are satisfied. 

The output from the Solution Procedure phase is the Results Transfer phase. This step in 
the methodology transfers the equilibrium output from the solution procedure to the final 
representation appropriate for post-processing and the crack growth model. 

1.5 Crack Growth Model 

The final step, and the essence of the crack growth methodology, is the Crack Growth 
Model. Figure 1.2 focuses in more detail on the Crack Growth Model. The input to the 
model is R/, the final representation, based on the current equilibrium state. The first 
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step of the Crack Growth Model is to apply the fracture and direction criteria to deter- 
mine the crack growth direction. In the case of linear elastic fracture mechanics the di- 
rection criteria may be dependent upon the fracture criterion (Erdogan and Sih, 1963), 
while in the case of elastic-plastic fracture mechanics the direction criterion is more 
likely to be dependant upon the stress or strain state local to the crack tip, and may or 
may not be tied to the fracture criterion in a theoretical sense. 

The Crack Extension Algorithms grow the crack using the predicted direction and user 
specified crack growth increment. This step of the methodology first deletes the mesh in 
the region local to a crack tip, then extends the crack. 

The Remeshing Algorithms place a rosette of elements around the crack tip, then rebuild 
the mesh around the crack tip. The crack extension and remeshing algorithms are essen- 
tially the same as for the original LEFM implementation, but modifications transparently 
manage elastic-plastic related data. 

The fourth step of the Crack Growth Model is essential for elastic-plastic crack growth. 
The Mapping Algorithms are responsible for maintaining the plastic history and dis- 
placement solution across the mesh changes created by the crack extension and remesh- 
ing algorithms. The plastic history includes primarily nodal displacements and gauss 
point stresses, but also includes the stress work, radius of the yield surface, and flags that 
indicate yielding state. 

The final step of the Crack Growth Model is to return the model to equilibrium, resulting 
in a new initial representation model, R k , at the new geometric configuration. Conceptu- 
ally, this step takes place as part of the crack growth model. For efficiency reasons, dur- 
ing automatic propagation equilibrium is re-established as the initial step of the solution 
procedure. 

1.6 Objectives and Scope 

The primary objective of this work is to develop a finite element model capable of mod- 
eling quasi-static, elastic-plastic, mixed-mode crack growth. A secondary objective is to 
implement the model in a relatively general way and in relatively user-friendly way, such 
that researchers and practicing engineers can apply the model to realistic problems. 
These objectives were met by extending an existing fracture simulation framework to in- 
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elude the necessary algorithms to model elastic-plastic fracture. The remainder of this 
work discusses the implementation of these algorithms. 

This work is neither a treatise on fracture mechanics, nor a treatise on finite elements. 
Rather, it addresses the issues necessary to model curvilinear elastic-plastic crack growth: 
a means of representing a crack without knowing a priori where the crack will go. In the 
process, discussions of supporting fracture mechanics and finite elements principles are 
included to provide context for modeling assumptions necessary to model the temporal, 
evolutionary nature of the problem. 

The Crack Growth Methodology illustrated in Figure 1 .2 essentially serves as the outline 
for the remainder of this work. Chapter 2 discusses various possible fracture criteria and 
direction criteria. Limitations of criteria in the literature are addressed, with the primary 
focus remaining on criteria that can be implemented on a framework that allows multiple 
interacting cracks in a single structure. Chapter 3 discusses crack extension algorithms 
for elastic-plastic analysis. The original algorithms were modified to allow necessary 
underlying data structures to manage elastic-plastic data. In addition, new algorithms 
were introduced to allow elastic-plastic crack growth along an edge of symmetry to allow 
symmetry when the crack lies on an edge of symmetry in the specimen. Chapter 4 dis- 
cusses the mapping algorithms. These algorithms preserve the plastic history as the mesh 
changes with crack growth. Chapter 5 presents analysis results for a number of Mode-I 
symmetry crack growth analyses and compares the results with experimental results from 
the literature. Chapter 6 presents analysis results for both Mode-I and mixed-mode crack 
growth using crack propagation with remeshing and mapping. The Mode-I results are 
compared with previous symmetry analysis results as well as experimental results. The 
mixed-mode analysis results are compared only with experimental results from the lit- 
erature. Chapter 7 summarizes this work and makes recommendations for future work. 
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Figure 1.1: Conceptual Crack Growth Methodology: how the crack growth model fits within the 
analysis architecture (after Potyondy, 1993). 



Figure 1 .2: Crack Growth Model: sequence of steps to grow a crack. 
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Chapter 2 


Fracture and Direction Criteria 


Two essential components of a software framework to predict mixed-mode tearing are a 
fracture criterion and a direction criterion. This chapter discusses a number of alterna- 
tives for both fracture and direction criteria, from both an historical perspective and a 
materials perspective, to provide insight into the decisions made for the criteria imple- 
mented. 


2.1 Fracture Criteria 

Early work in fracture mechanics focused on both the energy release required to create 
new crack surface area (Inglis, 1913), as well as the stresses in the body assuming iso- 
tropic linear elastic material behavior (Westergaard, 1939; Sneddon, 1946; Irwin, 1957; 
Williams, (1957). A typical presentation of the result is that the stress in any linear elas- 
tic cracked body can be described by the following equation in polar coordinates with the 
origin at the crack tip 


a u = 


■VT. 




m = 0 


( 1 ) 


where oy is the stress tensor, r and 9 define a point in polar coordinates relative to the 
crack tip, k is a constant, and fy and gy are functions of 9. The higher order terms depend 
on geometry, but the first term is proportional to 1 / 4r ■ As r approaches zero, the first 
term approaches infinity, but the higher order terms remain finite or approach zero. The 
result is that the first term dominates as r approaches zero. The stress intensity factor, K, 
is a convenient measure of the amplitude of the stress near the crack tip, where 
K = k-Jlx . Figure 2.1 shows the three modes of loading that can be applied to a crack. 
In Mode-I loading the load is applied normal to the crack plane; in Mode-II the loading 
corresponds to in-plane shear, which slides the crack faces relative to one another in the 
crack plane; and for Mode-Ill the loading refers to out-of-plane shear. The components 
of the three modes can be superimposed for the total solution of the problem. For a plate 
with a horizontal crack and a uniform far field stress perpendicular to the crack only 
Mode I loading is present. In this case equation (1) reduces to 
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( 2 ) 


K, 

w ^2 nr 

Figure 2.2 is a schematic plot of the stress normal to the crack plane due to the truncated 
form of Equation (1) compared to the stress due to the full expansion of the stress. The 
figure shows that for small values of r near the tip, the solutions match. This region 
where the solutions match is typically called the singularity dominated zone or the K 
dominated zone. In this region the stress-intensity factor completely defines the crack tip 
conditions. According to Anderson (1995), this single parameter description of crack tip 
conditions turns out to be one of the most important concepts in fracture mechanics. 

For linear elastic fracture mechanics the stress intensity factor is well established as a 
fracture parameter; however, defining a single parameter to characterize elastic-plastic 
fracture mechanics has proven more difficult. For elastic-plastic problems two widely 
recognized and accepted approaches to characterizing fracture are the J-path integral and 
the crack-tip opening angle (CTOA), or equivalently, crack-tip opening displacement 
(CTOD). These criteria are representative of two classes of criteria: resistance curve 
criteria and constant critical parameter criteria. The remainder of this section first dis- 
cusses the ./-Integral as a fracture criterion because of its recent intense scrutiny and 
popularity. Then other resistance curve criteria are discussed. Finally, the crack tip 
opening displacement is discussed. 

2.1.1 J-Path Integral 

The J-integral was introduced by Rice (1968) and is defined as a path-independent line 
integral enclosing the crack tip: 



where w is the strain energy density, T, are components of the traction vector along the 
boundary of the contour, u, are the displacement vector components, and ds is a length 
increment along the contour F Rice defined the J-integral based on the deformation the- 
ory of plasticity, so the integral is strictly only defined for linear or nonlinear elastic sol- 
ids. Rice showed that J is the energy release rate in a nonlinear elastic body that contains 
a crack. For stationary cracks and monotonically increasing loads, the deformation the- 
ory of plasticity will accurately model the plasticity behavior. However, for growing 
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cracks a region of elastic unloading exists in the crack tip wake and nonproportional 
loading occurs in front of the crack tip. 

Hutchinson (1968) and Rice and Rosengren (1968) (HRR) independently showed that J 
characterizes the crack tip conditions and that when the stresses and strains vary as 1/r 
near the crack tip J remains path independent. But the HRR singularity, as the 1/r singu- 
larity at the crack tip has come to be known, displays an anomaly; that is, it predicts infi- 
nite stresses as r ->0. 


Experimental results show that for many high fracture toughness materials, plastic de- 
formation causes an initially sharp crack to show blunting (Wells, 1961). McMeeking 
and Parks (1979) performed finite element analyses of crack tip configurations that in- 
cluded large strain theory and finite geometry changes. They showed that the HRR sin- 
gularity, which neither considers the effect of the blunting of the crack tip on the stress 
fields, nor takes into account the large strains present near the crack tip, is invalid in a 
small region near the crack tip. Figure 2.3 shows some results of McMeeking and Parks 
(1979). Note that the axes are nondimensionalized and are invariant as long as the plastic 
zone is small compared to specimen dimensions. The graph shows that the results of the 
large strain analysis peak for xaJJ at approximately unity, and decreases as x approaches 
0. The HRR singularity is invalid within this region near the tip, where the stresses are 
influenced by large strains and crack blunting. 


Hutchinson and Paris (1979) showed that for fully yielded conditions, ./-controlled crack 
growth is valid when: 


b dJ 


J 


Ic 


—— = co» 
da 


1 , 


(4) 


where b is the net section ligament of the specimen. Further, for J to be a useful fracture 
parameter, Anderson (1995) notes that steady state crack growth is necessary, and that 
even when J-controlled crack growth is obtained, steady state values of J are rarely 
reached. The conclusion is that under many, if not most, situations, the ./-resistance curve 
is not independent of specimen configuration and so loses both generality and effective- 
ness as a fracture parameter. 
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2.1.2 Other Resistance Curve Criteria 

Despite the theoretical limitations of the ./-Integral as a fracture criterion, it continues to 
see wide spread use, primarily as a resistance curve approach to predicting fracture. 
However, a number of other fracture criteria have been reported in the literature and have 
been used with varying degrees of success. 

Dawicke and Newman (1997) evaluated seven different fracture criteria for Mode-I 
fracture. They grouped the criteria into two categories: (1) constant critical parameters 
and (2) resistance curve parameters. The constant critical parameter approach assumes 
that an event such as fracture, growth, or link-up (depending on the criterion), will occur 
when the critical value is satisfied. The resistance curve approach assumes a unique re- 
lationship between the resistance curve fracture parameter and the amount of crack 
growth. The criteria evaluated were: 

• Crack-tip opening angle (CTOA) (constant critical parameter) 

• Plastic zone link-up model (constant critical parameter) 

• Stress-intensity factor resistance curve (ATr) 

• Effective Stress-intensity factor resistance curve (Krm) 

• ./-integral resistance curve (Jr) 

• Crack opening displacement resistance curve (<^0 

• T*-integral resistance curve 


While the elastic-plastic resistance curve criteria could predict the Mode-I fracture be- 
havior, they are limited in that adequate resistance curve data must be available, other- 
wise long crack simulations will require curve extrapolation. In addition, the resistance 
curves were not widely transferable between specimen types. The CTOA constant criti- 
cal parameter criterion was able to accurately predict the behavior of all of the experi- 
ments, including single crack and multiple crack panel specimens. 

Other criteria have been evaluated in the literature for mixed mode fracture. Two of the 
most recent are included here. Lee et al. (1996) used a linear relationship between plastic 
energy and crack size to predict fracture. However, it’s not clear how to apply the crite- 
rion in the case of multiple cracks in the same body, especially if the crack tips are ap- 
proaching for link-up. Similarly, Kuang and Chen (1997) used the plastic strain energy 
and showed that the energy required for initiation is relatively independent of mode 
mixity. Energy criteria such as these show promise as mixed-mode fracture parameters. 
But in addition to the limitations inherent as resistance curves, they also have limitations 
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in multiple cracking scenarios. As long as the plastic zones for each of multiple tips are 
distinct, each has its own unique plastic energy fracture parameter. However, after the 
plastic zones for two tips link-up it is unclear whether the plastic energy for the two can 
be considered separately. Traditionally, a plastic zone link-up criterion would be used as 
a secondary criterion. Dalle Donne and Doker (1997) present <5 r mixed-mode experi- 
mental results which show that crack opening displacement resistance curves do not cor- 
respond with Mode-I results from either C(T) or M(T) specimens, so Mode-I test results 
may not apply for mixed-mode situations. 

To summarize the discussion on the resistance curve approach to characterizing fracture, 
there are resistance curve parameters that will work for characterizing mixed-mode frac- 
ture for a single crack under specific conditions. However, all resistance curve parame- 
ters have a limitation on crack size, and it is not clear how to apply global energy ap- 
proaches to multiple turning cracks. Resistance curves will continue to be used because 
of their ease of use, but they lack generality in that they are not widely transferable be- 
tween specimen configurations or specimen types. 


2.1.3 Crack Tip Opening Displacement 

The crack tip opening displacement (CTOD), or its equivalent, the crack tip opening an- 
gle (CTO A) was identified as a possible measure of fracture toughness by Wells (1961) 
when he noticed, during testing of high toughness materials, that the degree of crack 
blunting increased in proportion to the toughness of the material. Typically CTOA and 
CTOD are defined in terms of each other as: 


y/ = 2 tan -1 


2d 


(5) 


where y/ > y/ c is the fracture criterion for CTOA and S > S c is the fracture criterion for 
CTOD, both defined a distance d behind the crack tip. 


Dawicke and Sutton (1993) provide a summary of significant CTOA/CTOD work up 
through 1993 for 2024-T3 aluminum. Experimental data (Demofonti and Rizzi, 1991; 
Newman et al., 1992; Dawicke and Sutton, 1993) indicate that surface measured values 
of CTOA decrease to a constant value after a small amount of crack growth (approxi- 
mately l-2x the sheet thickness). Figure 2.4 shows experimental values of CTOA meas- 
ured at various crack lengths using optical microscopy and digital image correlation 
methods (Dawicke and Sutton, 1993). The data clearly show the decrease in CTOA for 
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short crack lengths in the range of l-2x the thickness of the material. After the initial 
transient region, CTOA is essentially constant for the remaining crack growth. Recent 
three-dimensional residual strength analyses (Dawicke and Newman, 1998) have shown 
that for specimens of the same thickness, a single value of the critical fracture parameter 
can predict residual strength over a wide range of specimen sizes and configurations for 
both the M(T) and C(T) specimens. 

Computational studies of local crack tip fields for stationary cracks in thin sheets (Horn 
and McMeeking, 1980; Newman et al., 1993a) have shown that substantial constraint de- 
velops in thin materials prior to crack growth, and that the constraint at mid-thickness is 
about twice the constraint that occurs at the surface. Other computational studies of local 
crack tip fields for moving cracks have shown that the local constraint stabilized after a 
small amount of crack growth (Sommer and Aurich, 1991). This is in accordance with 
the data in Figure 2.4. 

Dawicke and Sutton (1993) have shown experimental results that correlate with the com- 
putational studies of constraint. They report that the tearing process is 3D, both during 
the initial transient, and after CTOA has stabilized. The crack front profile can be curved 
and shows tunneling effects (even for specimens as thin as 2.3 mm), or the fracture sur- 
face can be slanted through the thickness, depending on the crack orientation with respect 
to the material rolling direction. Both the results of Newman et al. (1993a) and Dawicke 
and Sutton (1993) show that displacements on the crack front near the free surface are 
higher than in the interior of the specimen. Both the experimental data and the analyses 
discussed here support the idea that even for thin-sheet materials, significant constraint 
develops around the crack tip. In addition, both indicate that the constraint stabilizes af- 
ter a small amount of crack growth. 

However, 2D simulations of tearing using CTOA as the fracture criterion continue to 
show good correlation with experimental data. Figure 2.5 is an example of load vs. crack 
extension from Dawicke and Sutton (1993) for 2024-T3 aluminum. The figure compares 
two computer simulations with the experimental data - the first is for a plane stress 
analysis and the second is for a plane stress analysis with a plane strain core in the finite 
element mesh to account for the higher constraint at the crack tip. Figure 2.6 shows a 
schematic of the plane strain core concept (Newman et al., 1988). Dawicke and Sutton 
(1993) tuned the height of the plane strain core to obtain numerical results that match ex- 
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perimental data. After accounting for the constraint at the crack tip, excellent agreement 
was obtained (see Figure 2.5). 

Amstutz et al. (1995a, 1995b) used the crack tip opening displacement to characterize the 
mixed-mode fracture behavior of 2024-T3 and considered the opening and sliding dis- 
placements separately. They found that CTOD measured 1 mm behind the crack tip, af- 
ter an initial transient region, approached a constant value for all modes of loading, indi- 
cating that CTOD can be used as a mixed-mode fracture criterion for 2024-T3. They also 
showed that the critical value of CTOD is a function of the material rolling direction, 
varying by about 25 percent for a crack growing against the rolling direction, compared 
to growing with the rolling direction. 

In summary two different classes of failure criteria have been considered: resistance 
curves and constant critical parameters. Resistance curves have an inherent limitation on 
the maximum crack growth without extrapolation. In addition, is not clear how to apply 
most resistance curves for Mixed-Mode cracking. In one case, the resistance curve was 
not the same for Mixed-Mode as for Mode-I. When multiple cracks are involved, crack 
tip link-up may be problematic, since as two tips approach, energy based methods may 
prove difficult to evaluate. In contrast, the CTOD constant critical parameter has been 
shown experimentally to be constant over wide range of mode-mixity. Since CTOD is 
measured local to the crack tip, multiple crack interaction effects are expected to be neg- 
ligible (in Mode I, link-up was accurately predicted). Finally, CTOD is relatively simple 
to implement compared to domain integral methods. 

2.2 Direction Criteria 

For linear elastic fracture mechanics, the fracture criterion can be intimately tied to the 
direction criterion. The conceptually simplest criterion is the cre-max criterion, where the 
crack propagates in a direction perpendicular to the direction of maximum tensile stress 
(Erdogan and Sih, 1963). Stability and direction are determined based on well-known 
hypotheses of brittle fracture in the K\/K u plane. Linear elastic fracture mechanics pre- 
dicts reasonably well the behavior of brittle materials, partly because these materials 
show a relatively consistent behavior from pure Mode-I loading to pure Mode-II loading 
(Erdogan and Sih, 1963). Under pure Mode-I the crack goes straight ahead. As increas- 
ing mode-II is introduced for an initial crack geometry, kinking at the crack tip is intro- 
duced and continually grows with increasing Mode-II. For example, Figure 2.7 shows 
experimental results for the initial crack extension angle in polymethylmethacrylate 
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(PMMA). Also included is the predicted angle based on the maximum tensile stress cri- 
terion. The theory matches well with the material behavior over the full range of mode- 
mixity. 

In contrast, ductile materials do not always show this type of behavior because ductile 
materials have complex crack tip behavior at failure that is not easily characterized. Re- 
sults by Maccagno and Knot (1992) on HY130 lightly tempered steel tested at room tem- 
perature showed a dramatic deviation from traditional LEFM behavior. Sharp cracks in 
Mode-I loading failed along planes of maximum shear, which leads to "zig-zag" fracture 
profiles that cannot be described by traditional LEFM fracture criteria. In effect, while 
remote loading was Mode-I with respect to the crack, the local failure was Mode-II. In- 
terestingly, HY130 steel in the identical condition but tested at -196° C, can be described 
by the maximum tensile stress criterion under mixed-mode loading. 

Recent studies by Amstutz et al. (1995b) on 2024-T3 aluminum and by Dalle Donne and 
Doker (1997) on 2024-T3 aluminum and on 550 structural steel show that for these mate- 
rials there is a sharp transition from Mode-I type fracture to Mode-II type fracture. In 
Mode-II, a “shear crack” forms as a result of shear of localization at the crack tip. Dalle 
Donne and Doker (1997) described this behavior, based on experimental and numerical 
observations, in terms of two competing mechanisms acting at the crack tip for a ductile 
material under mixed-mode loading. Figure 2.8 depicts these two mechanisms schemati- 
cally for an initially smooth notch tip. Under mixed-mode loading one side of the notch 
blunts from tensile stresses, and the other side sharpens from shear strains. Mode-I fail- 
ure occurs when the tensile stresses in the blunt region dominate. Mode-II failure occurs 
when shear strains localize along the sharpened side of the notch. The competition be- 
tween the mechanisms is resolved in the way that the material micro-structure responds 
to the applied stress and strain fields. The PMMA always fails due to tensile stresses at 
the crack tip. The HY130 steel always fails due to shear localization at the crack tip. 
Other materials, such as the 2024-T3 and 550 structural steel show a transition between 
tensile dominated and shear dominated failure. 

Ghosal and Narasimhan (1994) provided finite element results that help explain these two 
failure mechanisms under mixed-mode loading. Their plane strain model assumed large 
strain, but small-scale yielding. A Gurson type constitutive model was used to account 
for micro-void nucleation, growth and coalescence. Their results also showed that a por- 
tion of the notch surface blunts and the remaining portion sharpens, and that the plastic 
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shear strain localizes in the sharpened portion of the notch. The micro-void volume frac- 
tion was always nearer the blunted part of the deformed notch and independent of mode- 
mixity. The computer model did not include a means to explicitly represent material 
failure in the shear band (micro-voids cannot form in the shear band because of the lack 
of necessary stress triaxiality), but shear failure was simulated using an assumed maxi- 
mum shear strain criterion. They considered a range of loading angles from pure Mode-I 
to pure Mode-II, using T= tan' 1 K\IK\\ to characterize the loading angle at the crack tip. 
Their results showed that for 30° < W < 90°, micro-void coalescence at the blunted por- 
tion of the notch is the dominant failure mechanism, and that for 0° < < 15° the mate- 

rial failed due to shear localization at the sharpened portion of the notch tip. For 15° < V 7 
< 30° both of the mechanisms were active, and cause failure simultaneously. 

These ranges of failure modes match reasonably well with the experimental results of 
Amstutz (1995b), and numerical results of Sutton et al. (1997). Posing their results in 
terms of mode-mixity, 7 / = tan' 1 K\IK\\, where the stress-intensity factors are taken from 
Sutton et al. (1997) for each loading angle, micro-void coalescence dominates for 40° < 
*P< 90°, and shear localization dominates for 0° < *F< 40°. 

Several different criteria have been discussed in the literature for predicting the direction 
of crack growth under elastic-plastic conditions. Mahanty and Maiti (1989) showed that 
the maximum tensile stress criterion will predict the crack growth direction for tensile 
dominated fracture. Maccagno and Knott (1992) indicate that it may be possible to ob- 
tain a relationship between the opening and sliding components of crack tip opening dis- 
placement at initiation to predict the transition between Mode-I and Mode-II growth. But 
to date, no such relationship exists. Kfouri and Brown (1995) and Kfouri (1996) present 
a direction criterion based on maximum energy release rates at the tip of a short kinked 
crack. They develop the theory for LEFM, but show that a transition from tensile to 
shear failure can exist dependent on material parameters. Lee et al. (1997) used the 
maximum tensile force at the crack tip node. 

Sutton et al. (1997) analyzed the Amstutz (1995a, 1995b) data using small strain, small 
scale yielding, but with two terms in the theoretical asymptotic stress-intensity factor ex- 
pansion. He also described this type of ductile material behavior in terms of void initia- 
tion and growth concepts for the aluminum alloy 2024-T3. Sutton found that the direc- 
tion of propagation could be described in terms of a ratio of the mean stress cr m , and the 
effective stress <r e . When the ratio is greater than critical, the propagation angle corre- 
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sponds to the direction where the ratio is maximum; if the ratio is less than critical, the 
propagation angle corresponds to the direction where the effective stress is maximum. 
The mean stress and effective stress are evaluated on a circle of small radius around the 
crack tip. When the ratio is large, a high level of stress triaxiality exists. This empirical 
approach is attractive, but implementation details would be problematic. Finding maxi- 
mum values of <r m and cr e require “circle-plot” evaluations of these stress functions at 
small radii around the crack tip. Unfortunately, the maximum values are located on wide 
plateaus where very small numerical errors can lead to erroneous results. 

Most recently, Chao and Liu (1998) considered a ratio of the maximum shear stress to 
maximum tensile stress, evaluated on a circle near the crack tip. They showed that a 
critical value of the ratio exists that predicts the transition from shear to tensile failure. If 
the ratio is less than critical, tensile failure exists. If the ratio is less than critical, shear 
failure exists. As with the Sutton et al. (1997) criterion, implementation may be prob- 
lematic since results evaluated on a circle near the crack tip may be on wide plateaus 
where very small numerical errors can lead to erroneous results. 

In summary, experimental evidence suggests that different materials have different domi- 
nant failure mechanisms. No clear theoretical basis as yet is available in the literature 
that explains the transition between tensile failure and shear failure. It is possible that an 
empirical transition based on experimental evidence, coupled with a dual failure mecha- 
nism, maximum tensile stress, and maximum shear stress could predict the direction, but 
such a criterion was not investigated. A wide range of materials have most of their be- 
havior characterized by the maximum tensile stress criterion. In addition it is somewhat 
unusual (though not impossible) for cracks to initiate and grow under Mode-I loading, 
then be subjected to an extreme overload in Mode-II. That is to say, most cracks start 
and grow dominated by Mode-I, with curving to minimize Mode-II. 

2.3 Summary 

In summary, the fracture criterion and direction criterion are inextricably tied to the crack 
tip stress and strain state, but no convenient, unified theory exists that defines the fracture 
behavior for elastic plastic materials. The above referenced recent research indicates that 
a number of different fracture criteria and direction criteria can be used under certain cir- 
cumstances. For the work presented here the fracture criterion is assumed to be the mag- 
nitude of the crack tip opening displacement vector, and the direction criterion is assumed 
to be maximum tensile stress at the crack tip. CTOD is largely unproven to be constant 
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under mixed-mode loading for materials other than 2024-T3 Aluminum, so its applica- 
bility to other materials is somewhat in question. It’s clear that the maximum tensile 
stress criterion will not correctly predict the fracture direction for all cases. However, the 
criterion is simple, requires no special materials parameters for implementation, and does 
perform well under a wide variety of circumstances for a wide variety of materials. More 
widely applicable criteria are an area for future research. 
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Figure 2.1: The three modes of loading that can be applied to a crack (Anderson, 1995). 



Figure 2.2: Schematic of stress singularity zone (Anderson, 1995). 


19 



Figure 2.3: Large strain crack tip finite element results of McMeeking and Parks (1979). Blunting 
causes the stress to deviate from the HRR solution close to the crack tip (Anderson, 1995). 



Figure 2.4: Typical experimental values of crack tip opening angle measured for 76.2 mm wide, 
2.3 mm thick M(T) specimens, where B is the specimen thickness (Dawicke and Sutton, 1993). 
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Figure 2.5: Load vs. crack extension data shows the effect of the plane strain core for 2024-T3 
aluminum (Dawicke and Sutton, 1993). 



Figure 2.6: Schematic of the plane strain core. 
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Figure 2.7: Experimental results for the initial crack extension angle in PMMA (Maccagno and 
Knott, 1989). 
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Figure 2.8: Competing fracture mechanisms at a crack tip (Dalle Donne and Ddker, 1997). 
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Chapter 3 


Crack Propagation Algorithms 


For more than two decades special purpose finite element codes have been used to model 
crack growth in the presence of plasticity. This chapter describes many of these models 
and then outlines the implementation and capabilities of the algorithms implemented for 
this work. 

3.1 Background 

The earliest crack propagation algorithms were based on a nodal release technique, 
whereby the crack path is predefined in the finite element mesh (Anderson, 1973; de 
Koning, 1977; Newman, 1974). When the fracture criterion is satisfied or exceeded the 
forces imposed to hold two coincident nodes at a crack tip are released, thus creating new 
surface area for the crack and a new crack tip at the adjacent node. Typically, stiff spring 
elements are used to tie the coincident nodes together. This technique is still in use today 
in both two- and three-dimensional codes (Maiti and Mahanty, 1990; Roy et al., 1993; 
James, 1997; Dawicke and Newman, 1998). Nodal release is typically applied to Mode-I 
cracks, but can be applied to curving cracks with a predefined crack path. 

A number of other approaches to modeling crack propogation have also been published 
in the literature, and representatives of these are included here. Nakagaki et al. (1979) 
analyzed Mode-I crack growth by deforming the crack tip elements such that the crack tip 
shifts ahead and employed a mapping procedure to transfer state variables to the modified 
mesh however remeshing was not included. Deng (1990) studied crack tip fields for dy- 
namic fracture. His software assumed small-scale yielding at the crack tip, but included 
effects of the plasticity. An Eulerian formulation allowed the finite mesh to translate in 
steady state with the crack tip. The model included effects for only Mode-I loading. Lim 
(1992) modeled mixed-mode fracture using discrete cracks with automatic remeshing 
where there was plasticity in the crack wake, but the direction and fracture criterion were 
governed by LEFM. Mapping was used to transfer state variables for the remeshing re- 
gion. However, since remeshing occurred only in the LEFM region in front of the 
growing crack, only displacements were mapped. The technique by Bakuckas et al. 
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(1993) does not require a predefined crack path, but requires the crack path to follow ex- 
isting element boundaries in the mesh using a node splitting technique. 

More recently, Marusich and Ortiz (1995), implemented a dynamic Lagrangian finite 
element program to predict fractures of high-speed machining processes. The model in- 
cludes both a shear localization mechanism to initiate fracture and a traditional stress in- 
tensity based fracture criterion to propagate cracks. A continually updated mesh is 
maintained by automatic remeshing algorithms. 

Camacho and Ortiz (1996) report on a similar dynamic Lagrangian finite element pro- 
gram applied to modeling impact damage in brittle materials. The fracture criteria are 
based on maximum shear and maximum tensile stress values. Cracks are restricted to 
propagating along existing element boundaries and use shear and tensile strain softening 
behavior. 

Lee et al. (1997) recently published work that is the most directly related to this work. 
The authors modeled non-self-similar crack growth using finite elements and small strain 
plasticity. The fracture criterion is based on a global plastic work criterion, whereby the 
crack propagates according to a linear relationship between plastic work and crack exten- 
sion. The direction criterion has the crack propagate perpendicular to the maximum ten- 
sile force at the crack tip node. The actual crack propagation algorithms used were not 
clearly specified. It appears that at each propagation step an external program is used to 
determine propagation direction. Then a separate external program somehow updates the 
mesh in the new configuration. Finally the analysis is resumed using load relaxation to 
relax crack face tractions before external loading resumes. 

While this most recent approach is similar to the approach used here, some differences 
exist: 

• In the current work a local fracture criterion is used, while Lee et al. (1997) use a 
global energy fracture criterion. The local approach makes possible multiple cracks, 
each following a different curving path in the model. The energy approach is untested 
for curving cracks, and is unlikely to work for multiple cracks, since only the global 
energy is considered, not the energy at each tip (currently the T* integral shows prom- 
ise as a local crack tip energy-based criterion (Dawicke and Newman, 1998)). 
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• The current implementation for this work can include multiple cracks in multiple lay- 
ers and has fully automatic crack propagation without user intervention. 

3.2 Combined Remeshing/Nodal Release Algorithm 

No comprehensive algorithms are referenced in the literature that can handle multiple 
curving and interacting cracks in the presence of plasticity (non-small scale yielding con- 
ditions, net section yielding, etc.). The approach adopted for solving this class of elastic- 
plastic fracture problems is to extend the LEFM approach of Wawrzynek and Ingraffea 
(1987) to include the effects of plasticity. Their approach to fracture simulation uses an 
explicit representation for each crack. Crack growth is modeled using remeshing algo- 
rithms that update the mesh only locally to the propagating crack. For LEFM analyses, 
no history variables are necessary, and a new linear analysis is run at each propagation 
step. For elastic-plastic fracture mechanics (EPFM), the plastic history is important both 
for both the plasticity in front of the crack tip and the residual stresses in the wake; this 
history must be accurately preserved even as the mesh is continually changing with crack 
growth. State variable mapping is a typical approach to preserving the history as the 
mesh changes. Critical element and node quantities are saved before the remeshing proc- 
ess starts. After remeshing is complete, these quantities are mapped onto the new mesh 
elements and nodes using the shape functions of both the old and new mesh. State vari- 
able mapping has been used extensively for other applications such as metal forming and 
multi-grid finite difference analyses. The specifics of state variable mapping will be 
elaborated upon in a later section. 

The local remeshing and state variable mapping approach to EPFM has been combined 
with the nodal release technique and implemented in the FRANC2D/L program to pro- 
vide a powerful and flexible tool to model elastic plastic crack growth. The nodal release 
technique is implemented such that crack propagation along lines of symmetry can be 
modeled for optimum mesh efficiency, and the local remeshing algorithms have been 
augmented with state variable mapping algorithms to allow fully automatic elastic plastic 
tearing without user intervention. In addition, the two algorithms can be combined to al- 
low the nodal release technique along a predefined curving crack path. As an added 
benefit for flexibility, the local remeshing technique can also be used in a step-by-step 
mode, where the crack path is user specified. The following sections briefly describe each 
of the four variations on how elastic-plastic fracture can be modeled before continuing on 
with a more detailed description of the algorithms. 


25 



3.2.1 Nodal Release 

Symmetry Edge Many analyses of laboratory specimens will have natural lines of sym- 
metry that can be exploited during the analysis to reduce the mesh size by half or more. 
This option allows crack propagation along the edge of symmetry that follows a hori- 
zontal axis. Figure 3.1 is a schematic showing how lines of symmetry reduce the model 
size of a middle crack tension specimen to 14 of the full specimen size. Symmetry 
boundary conditions along the vertical centerline of the specimen model the side to side 
symmetry. Since there is no cracking along the vertical centerline, these boundary con- 
ditions will not change during the analysis. Symmetry boundary conditions also exist 
along the crack face. The part of the crack face that is fully closed is held by fixity con- 
ditions that can be efficiently released as crack propagation becomes necessary. This 
form of crack propagation is commonly called the nodal release algorithm. In the current 
implementation, multiple symmetry cracks can exist and need not lie on the X axis, but 
they must all lie on a horizontal edge of the model. 

Predefined Path Under some circumstances, users may know the crack path in advance 
(e.g. from experiments) and want to force the analysis to use that path in the most effi- 
cient manner. For this case, the nodal release algorithm can be applied for a predefined 
curving crack path. Figure 3.2 is a schematic showing nodal release along a predefined 
path. The crack has coincident node pairs on each side of the crack from the crack mouth 
to the crack tip. Coincident nodes are efficiently tied together in a manner similar to the 
symmetry edge method. 

3.2.2 Local Remeshing and Mapping 

User Specified Path Under some circumstances, users may want to monitor the crack tip 
state manually at each propagation step to extract response data to implement a custom 
direction criterion. In this case, the user can explicitly specify the crack path. State vari- 
ables are mapped at each remeshing step to insure a correct nonlinear history for the 
yielded material. 

Fully Automatic Propagation This option provides the user with automatic propagation. 
At each propagation step the crack growth direction is determined from a direction crite- 
rion, and the step size is user specified. Automatic remeshing algorithms extend the 
crack and update the mesh at the new crack-tip location. State variables are mapped at 
each remeshing step. 
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3.2.3 Algorithm Description 

The capabilities desired - modeling elastic-plastic fracture with the level of flexibility 
described previously - led to two basic criteria for the implementation: relative simplicity 
for the user to specify a problem, and flexibility and simplicity in the internal database 
structure of the underlying software. These criteria were met by designing a unified data 
structure for the implementation and by establishing automated procedures to maintain 
the data structures in a manner transparent to the user. 

However, the resulting implementation is still rather complex. All propagation algo- 
rithms revolve around the nodal release algorithm, since that implementation is used to 
monitor CTOD and other crack tip fracture criteria during the analysis. In addition, the 
nodal release algorithm is embedded within the nonlinear Newton solution scheme, since 
an efficient implementation of nodal release is intimately tied to the Newton scheme. 
However the remeshing algorithms significantly modify the mesh, and require that the 
system tangent stiffness matrix be rebuilt from scratch. Listing 3.1 shows the conceptual 
model for the automatic tearing algorithm without regard to how these capabilities are 
implemented. The steps are numbered only for clarity as part of this discussion. 


> '■ 

2. 

-> 3. 

4. 

5. 


6 . 

7. 



8 . 


Initialize 

Begin analysis using current state 
Increment load/displacement BCs 
Find equilibrium 

Check before tearing. If detected: 

If unzipping, untie crack face nodes 
If re meshing, prepare to return for remesh 
Until total load 
Petformed remesh if require 
7a. None required-done with tearing analysis 
7 b. Save state for full mesh 
7c. For each crack tip: 

Delete locally 
Extend crack 
Add new mesh 
7d. End loop 

7e. Map saved state to new mesh 
7f. Return to equilibrium 
Until tearing analysis done 


Listing 3.1: Steps of the automatic tearing algorithm. 
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Step 1 of the algorithm is initialization. This includes not only general initialization 
steps such as boundary conditions and analysis control data for the non-linear solver, but 
also initialization specific to the tearing algorithm. For instance, internal cracks that tear 
with remeshing require different user input for initialization than edge symmetry cracks 
that tear with nodal release. 

Step 2 is to begin an analysis using the current state. This is an important conceptual 
consideration. The current state may be and initial non-loaded model or a restart state. A 
restart state may be the result of restart from a file, the result of a user request to view re- 
sults, or the result of mapping algorithms after a request to grow a crack. From the per- 
spective of resuming an analysis, each of these is considered equivalent. 

Step 3 is to increment the load and displacement boundary conditions (BC). The current 
algorithm for incrementing the BC’s is quite simple. A total user defined load increment 
is broken down into a linear distribution of sub-increments (also user defined). Equilib- 
rium is established at each sub-increment. If sub-increments are too large, either the 
plasticity algorithm will fail or the Newton scheme will fail. In either case an error mes- 
sage indicates the nature of the problem. Automatic load stepping schemes are available 
in the literature (c.f. Crisfield, 1991) but are not implemented here. One exception is a 
simple, but effective predictor algorithm to track the fracture criterion and predict the 
load step fraction that allows the fracture criterion to be met (within a tolerance). 

Step 4 is to find equilibrium at the current loading condition. This entails using the tan- 
gent stiffness matrix in the modified Newton iterative scheme. Appendix A provides a 
more detailed description of the finite element derivation and implementation. 

Step 5 is to check if any of the current crack tips have a fracture criterion within tolerance 
of critical, and if so, take the appropriate action. In the case of nodal release cracks, the 
release is handled locally so that the analysis can continue uninterrupted. In the case of a 
crack tip that requires remeshing, the analysis is interrupted, and automatic remeshing is 
invoked (see Step 7). 

Step 6 is presented as looping until total load has been exhausted. This is true for simple 
plasticity analyses without cracking, but for tearing the termination criterion is the num- 
ber of crack propagation steps. 
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Step 7 is to perform remeshing if required. This step is executed after exiting the loop 
over steps 3-6, where the exit condition could be any one of several conditions, including 
an error condition or a request that remeshing take place. 

Step 7b is to save the state of the analysis for the full mesh. This step is part of the state 
variable mapping process. 

Step 7c is the remeshing process, where for each crack tip the mesh is deleted locally to 
the crack tip, the tip is extended, and a new mesh is created and added around the tip. 

Step 7e is to transfer the state saved in step 7b onto the new mesh using the state variable 
mapping routines. 

Step 7/ returns the newly mapped state to equilibrium. The state was mapped onto a new 
geometry that has new crack surface area. Consequently, the compliance of the structure 
will have changed. In addition, mapping errors are introduced as a result of approxima- 
tions in the mapping process. Both of these effects are somewhat accounted for by re- 
turning the state to equilibrium before resuming the analysis. (Mapping errors are dis- 
cussed in Chapter 4.) 

Step 8 simply indicates that the analysis will continue until the tearing is complete. As 
indicated in Step 6, the user specified maximum number of propagation steps controls 
analysis termination. 

3.2.4 Nodal Release Implementation 

The nodal release mechanism uses stiff springs to efficiently model crack growth for 
cases where the crack is growing along a symmetry edge or along a predefined path. 
Figure 3.1 is a schematic that shows the conceptual model for the nodal release mecha- 
nism. Stiff springs provide forces to hold coincident nodes together (or to provide zero 
displacement, in the case of symmetry edge). When the crack tip fracture criterion is 
met, the stiff springs at the tip are removed from the problem and corresponding crack 
face forces can be relaxed, while holding all other boundary conditions constant. 

Since the stiff springs are effectively part of the mesh, in the most general implementa- 
tion the full stiffness would need to be re-formed and re-factored. However, an imple- 
mentation by Newman (1974) can efficiently remove the springs from the stiffness matrix 
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in its inverted form without re-factoring. The algorithm is based on a modified Cholesky 
algorithm and greatly enhances the efficiency of the nodal release algorithm. 

The value of the spring stiffness constant is an important consideration. The effect of the 
stiff springs is to enforce a penalty function for the equations of the global system corre- 
sponding to the nodal displacements. The ideal effect would be to add the stiff spring 
before propagation and remove the stiff spring after propagation without effect on the 
global stiffness matrix. In practice, the value of the stiff spring needs to be large enough 
that the stiffness due to all other elements framing into the node is insignificant compared 
to the stiff spring value. However, when the stiff spring is subtracted out of the inverted 
stiffness matrix, the original global stiffness coefficients must be recovered. Since all 
operations are performed in double precision, there are about 15 significant digits avail- 
able. For example, if a coefficient has a value of 1.0E4 and the stiff spring with a value 
of 1.0E10 is an added, then subtracted from the coefficient, about six significant digits 
will be lost and about 9 accurate significant digits will be retained. Ideally the springs 
stiffness would be determined as a function of the global stiffness matrix. In the current 
implementation a constant K s —l -0E4 is used. 

3.3 Crack Tip Constraint Model 

Stress triaxiality is a significant issue at the tip of a crack, even for thin sheet material. 
This stress triaxiality, or constraint, has received much attention, but still few models ex- 
ist to describe constraint in a two-dimensional analysis. Plane stress has no constraint, 
and plane strain allows the constraint to extend too far away from the crack tip. 

Newman et al. (1993a) describes constraint as a buildup of stresses around a crack front 
due to the resistance against in-plane and out-of-plane deformation. In-plane constraint is 
mainly associated with the closeness of the crack front to external boundaries, and out-of- 
plane constraint is mainly influenced by plate thickness. While in-plane constraint is an 
important fracture issue, from the perspective of residual strength in-plane constraint is 
less of immediate concern, since in many cases, only a small amount of crack growth 
may be necessary before a structure experiences maximum load. In addition, by the time 
the crack is being significantly influenced by an external boundary, fracture may be con- 
trolled more by the stress-strain behavior of the material and net-section yielding in the 
remaining ligament. Out-of-plane constraint, on the other hand, can have significant ef- 
fects on the fracture behavior and residual strength of structures very early in the tearing 
history, before the maximum load is reached. 
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Anderson (1995) describes out-of-plane constraint relative to an uncracked plate sub- 
jected to in-plane loading. With no crack, the plate is in a state of plane stress. If the 
plate has a small through-crack perpendicular to the load somewhere far from all bounda- 
ries, the material far from the crack is still in a state of plane stress. The stresses are ele- 
vated near the crack tip, and these elevated stresses cause a Poisson effect along the crack 
front, causing the material to contract; the surrounding material resists the contraction, 
resulting in out-of-plane stresses along the crack front. That is, the crack is in a triaxial 
state of stress, while the material far from the crack is in a biaxial state of stress. 

Figure 3.3 shows a schematic of how the stresses and strains vary along the crack front. 
On the interior of the crack, the material sees effectively a state of plane strain. Where 
the crack intersects the free surface, the material is in a state of plane stress. A transition 
region exists near the free surface where the material is neither plane stress nor plane 
strain. As the material around a crack front yields, the stress state will vary along the 
crack front from plane stress at the free surface, to some amount of triaxial state at the 
mid-plane of the crack. Once the plastic zone size reaches approximately half of the ma- 
terial thickness, the stress state is predominantly plane stress in the yield region, except 
for very near the crack front. 

There are a number of ways to characterize constraint (Newman et al., 1993a). One sim- 
ple way is to look at the out-of-plane stress near the crack front. Figure 3.4 shows the 
normalized Z stress as a function of position in front of the crack at the mid-plane of the 
specimen. The figure shows that the Z stress, and thus constraint, decreases rapidly to 
zero in less than 1/2 of the thickness of the material. 

Newman et al. (1993a) used three-dimensional elastic-plastic finite element analysis to 
investigate the constraint along a stationary crack front. Newman's results showed that 
even under large-scale plastic yielding, local constraint developed in the interior regions 
near the crack even for the thinnest material investigated (B = 1 .25 mm). 

3.3.1 Constraint and Residual Strength 

Early attempts to simulate residual strength with plane stress or plane strain analyses 
were unsuccessful. Fracture parameter values “tuned” to match residual strength for a 
middle-crack tension specimen would not match analysis results for a compact tension 
specimen. The reason for the inconsistent analysis results was inability of either a fully 
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plane stress or fully plane strain model to capture both the plane stress state away from 
the crack and triaxial stress state near the crack. 

There are number of different ways to simulate out-of-plane crack tip constraint in a two- 
dimensional analysis. Shan et al. (1992) uses either three-dimensional analysis or ex- 
perimental results as a baseline, and correlates plane stress and plane strain results to get 
a two-dimensional approximation to the three-dimensional problem. The correlation is 
performed by running both plane stress and plane strain analyses, then using a linear 
combination of the two, assuming that a fraction of the thickness B s /B is in plane stress 
and that the remainder of the thickness is in plane strain (1-B S /B) (B s is the thickness of 
the plane stress section, and B is the total thickness). The approach produces reasonable 
results, but requires two analyses, with superposition after the fact. Shan et al. (1994) 
showed that his approach is somewhat independent of specimen size; however, his work 
was for straight cracks. It is unclear how the work would translate to mixed-mode be- 
havior and whether the approach is specimen independent. 

Newman et al. (1988) introduced a mixed state of stress, commonly called the plane 
strain core, to capture the constraint at the crack tip. Figure 3.5 is a schematic of the 
plane strain core concept illustrating Newman’s mixed state of stress. Plane stress ele- 
ments are used for the entire model, except along the crack face and in front of the crack 
tip, which are modeled with plane strain elements. The plane strain core has been used to 
accurately predict residual strength for C(T) and M(T) specimens with the nodal release 
algorithm (Newman et al., 1988; James, 1997). FRANC2D/L has this plane strain core 
capability for nodal release fracture simulations. 

Newman’s plane strain core did not evolve as the crack grew because the core was fixed 
in front of the crack line. The predefined core works well for problems for which the 
crack path is known in advance. However, when the crack path is not known a priori, the 
region of high constraint must move with the crack tip. In the three-dimensional analysis 
the constraint develops naturally near the crack as part of the analysis. 

3.3.2 Constraint and the Remeshing Model 

When the nodal release mechanism is used to simulate fracture the mesh does not change 
during crack growth. As result, a constant height plane strain core can be “tuned” to 
model the constraint at the crack tip, and the core can be predefined in the mesh. When 
the remeshing-mapping crack growth mechanism is used to simulate fracture the mesh 
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does change during crack growth. Since the crack path is not known in advance, the 
plane strain core cannot be predefined and must grow as the crack grows. Two sets of 
issues surround implementing a growing plane strain core: those related to representing 
the geometry of the core, and those related to the underlying continue mechanics. The 
implementation of a solution to these issues is not within the scope of this work. 

3.4 Simulating Multi-Site Damage 

One cracking phenomenon of particular concern to the aircraft industry is the residual 
strength of the structure in the presence of multi-site damage (MSD). MSD cracking in- 
volves two or more cracks that can ultimately interact to influence each other's behavior 
and the behavior of the overall structure. Two perspectives must be considered when 
simulating multi-site damage: the nodal release algorithm and the remeshing-mapping 
algorithm. 

3.4.1 Nodal Release MSD Model 

Multi-site damage in the context of the nodal release algorithm is conceptually simple. 
The nodal release approach ties nodes along an edge model. Figure 3.6 shows a sche- 
matic of an M(T) specimen with a main “lead” crack and MSD in front of the lead crack. 
All of the nodes along the bottom edge of the model were initially tied. Next, the lead 
crack was specified, which untied nodes from the middle of the specimen to the tip of the 
lead crack. Finally, nodes were untied at the specific locations of the MSD in front of the 
lead crack. 

As the analysis progresses, all of the tips will be checked for the fracture criterion, and 
nodes released as appropriate to allow crack growth along the crack face. The criterion 
for merging two cracks, such as two MSD cracks, or the lead crack and MSD crack, is 
simply that when only one node remains between two cracks, that node is released auto- 
matically, and the cracks are merged. Figure 3.6 also shows schematically the case 
where two MSD cracks are ready to merge. Only two nodes remained tied. When either 
of the two tips reaches the critical fracture criterion, the node at that tip will be released, 
leaving only one node attached between the two MSD cracks. When only a single mode 
remains that node is released automatically. 
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3.4.2 Remeshing-Mapping with Multiple Cracks 

When remeshing is used, the situation is somewhat different. Since the crack path is not 
predefined, the geometry ahead of each crack tip must be detected and handled robustly 
within the algorithms. Several possible situations are: 

Tip Approaching Tip. This fairly general situation has the previous M(T) example 
in Figure 3.6 as a special case. Figure 3.7 is a schematic of the situation for the 
remeshing analysis. In general the situation arises when two tips are close enough 
that the ligament between them fails. 

Tip Approaching Crack Face. In an analysis with curving cracks this situation is 
probably most likely. A classic example is the “football” formed by two fatigue 
cracks approaching, then repelling one another, then merging by intersecting the 
crack face. Figure 3.8 is a schematic of this type of crack coalescence. 

Tip Approaching Boundary. During some fracture scenarios, it may be possible 
for a crack to approach and intersect an outer boundary of the model before residual 
strength is reached. This condition is a concern because the current implementation 
does not allow multiple bodies (full fracture into two or more pieces). 

In the current implementation each of these situations is handled uniformly and simply, 
but not automatically. The algorithm does detect automatically when a tip approaches the 
boundary. At that point control is returned to the user. Manual remeshing procedures 
can be used to extend crack tips manually, or to delete elements to manually coalesce 
cracks. Manual mapping procedures can be used to retain the plasticity state and allow 
restart of an analysis. 

3.4.3 Load Relaxation 

When simulating MSD problems there are multiple crack tips. Any of these crack tips is 
free to grow during an increment in the boundary conditions. Under some circumstances 
more than one may grow during this increment. When the remaining ligament between 
two tips is small, the two tips can grow together, breaking the ligament. Under these cir- 
cumstances intermediate equilibrium states are recorded as the cracks are growing, pre- 
serving the plasticity history and extension state for each of the growing cracks. A load 
relaxation scheme is used, whereby the forces holding the crack closed before growth are 
slowly relaxed over ten steps, while holding the boundary conditions constant. If during 
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the load relaxation process load redistribution causes additional crack tips to meet the 
fracture criterion, the crack face loads at those tips are also relaxed over ten steps. The 
equilibrium iterations for relaxation continue until all crack face loads are expended and 
no more tips meet the fracture criterion. 

3.5 Summary 

This chapter described two fundamentally different approaches to elastic-plastic crack 
propagation: nodal release and automatic remeshing with mapping. The two approaches 
have been implemented in the FRANC2D/L code. The algorithms have been imple- 
mented to model elastic-plastic fracture under a number of different circumstances. The 
nodal release algorithm is specialized to work for either of two different situations: 
cracking along an edge symmetry or cracking along a predefined curving crack path. The 
automatic remeshing algorithm can be used either in a step-by-step manual mode or in 
fully automatic mode. In fully automatic mode the direction is predicted from maximum 
tensile stress at the crack tip. Crack tip constraint can be modeled during nodal release 
using the plane strain core. Specialized algorithms allow multi-site damage during nodal 
release. Manual remeshing and mapping can be used to merge MSD crack tips when 
automatic propagation is interrupted by crack tips in close proximity to one another. The 
algorithms implemented are flexible and powerful and can be applied in a wide variety of 
situations. 
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Figure 3.1 : Schematic showing lines of symmetry for a middle crack tension specimen. During 
elastic-plastic fracture, fixity maintain symmetry along the vertical edge, and special stiff springs 
maintain symmetry along the fracture surface. 


Path defined, 
Forces ahead of 
crack relaxed during 
propagation 



Figure 3.2: Schematic showing crack face forces for a nodal release fracture simulation with the 
crack path predefined. Crack face forces are maintained by special stiff springs. 
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Figure 3.3: Schematic showing the variation of stress and strain along a three-dimensional crack 
front (Anderson, 1995). 
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Figure 3.5: Plane strain core concept. 



Figure 3.6: Crack tip coalescence for a symmetry model. Multi-site damage configuration with 
tips ready to merge when criterion at either tip reaches critical. 
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Figure 3.7: Schematic of Tip Approaching Tip. Circles represent minimum radius for automatic 
propagation to continue. 



Figure 3.8: Schematic of Tip Approaching Boundary. Circles represent minimum radius for 
automatic propagation to continue. 
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Chapter 4 


State Variable Mapping 


The elastic-plastic constitutive equations are nonlinear and history dependent. For the 
finite element method, state variables are stored at the gauss points of each element to 
represent this history dependence. Typical state variables may include the stresses, hard- 
ening parameters, work variables, and yielding flags. After remeshing takes place, the 
state variables must be accurately transferred from the old mesh to the new mesh to retain 
the history information. In addition, nodal displacements must be transferred from the 
old node locations to the new node locations. 

4.1 Mesh Mapping 

This process of transferring state variables from one mesh to another is typically referred 
to as mesh mapping. Mesh mapping techniques are discussed in finite element literature 
for fracture simulations (Lim et al., 1992; Kato et al, 1993), as well as metal forming pro- 
cess simulations (Ponthot and Hogge, 1991; Sekhon and Chenot, 1993), welding proc- 
esses (Brown and Song, 1993), and in more general settings (Chavez, 1983). Four refer- 
ences that address issues concerned with mapping large strain response data are Ponthot 
and Hogge (1991), Kato et al. (1993), Yamada (1993), and Espinosa et al. (1996). Kato 
et al. (1993) note that accuracy of the displacement mapping is especially important 
“when using the total Lagrangian formulation, since the initial geometry must be cor- 
rectly recovered from the deformed geometry using the mapped displacements.” Since 
this work is only concerned with small strain and small geometry changes, these large 
strain issues will not be addressed further. 

The displacement mapping problem can be viewed as the process of transferring dis- 
placements from an old mesh on the current volume to a new mesh on the current vol- 
ume: 


U„(V)~> U„(V), 


( 1 ) 
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where V is the current volume in time, and U Q and U„ are old and new displacements cor- 
responding to the old and new mesh nodal coordinates, X„ and X„. And so, the mapping 
problem can be summarized as: 

Given U,„ X„ and X„, 

Find U„. 

By direct extension, the process can be applied to any nodal response variable. 

Several methods of transferring state variables between old and new meshes are dis- 
cussed in the literature. Lim (1992) summarizes the methods in three categories: 
smoothing techniques, triangulation, and inverse isoparametric mapping. The state vari- 
able mapping algorithm implemented for this work is the inverse isoparametric mapping 
method. 

4.2 Inverse Isoparametric Mapping 

Isoparametric elements use the same shape functions to interpolate the nodal coordinates 
and the nodal response data: 

x.^nm^-x:, < 2 > 

/ 

U„=J,N,(r,s)-Ul. 0) 

( 

where N, are the shape functions evaluated at the element natural coordinates (r,s), and X 0 
and U 0 are nodal values of coordinates and displacements, respectively. For the inverse 
mapping problem, given the coordinates of a node in the new mesh, find the natural co- 
ordinates in the old element, such that interpolation in the old element gives the coordi- 
nates of the new node. That is: Given X Q and X„, the nodal coordinates, find (r n ,s n ) such 
that 


x.^nm.'S.YK- < 4 > 

i 

Then ( r n ,s n ) can be used to interpolate the nodal response data from the old mesh into the 
new mesh: 
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( 5 ) 


/ 

Equation (4) is a nonlinear equation in ( r,s ) which cannot be solved directly. 

To solve Equation (4), Lim (1992) uses a predefined line technique, where a line is drawn 
across the element through the point where the natural coordinates are desired and to one 
of the nodes on the element. The natural coordinates are found iteratively by using a re- 
sulting set of simultaneous equations. The solution is straight forward for straight sided 
elements, but for curved sided elements many special cases exist. But when correctly 
implemented, Lim (1992) reports that the algorithm is robust and efficient. 

The current implementation uses the method of Kato et al. (1993), where the shape func- 
tion derivatives are used in a simple Newton method iterative algorithm. For straight- 
sided elements the algorithm converges rapidly and is quite accurate, and no special pro- 
visions are necessary for curved elements, should they become necessary in the future. 
The Newton method has the added benefit that it was already implemented and well 
tested in FRANC2D/L (used primarily for post-processing purposes in line plots of re- 
sponse data). 

4.3 Stress Mapping 

As previously mentioned, the inverse isoparametric mapping algorithm can be used to 
map any nodal response variable. Stresses and other element based response data must 
be mapped from gauss points in the old mesh, to gauss points in the new mesh. To ac- 
complish this several steps are necessary: 

1 . Gauss point global coordinates are evaluated and saved for gauss points in the new 
mesh using the gauss point natural coordinates. 

2. Gauss point response data from the old mesh are extrapolated to the element nodes 
using a least squares fit of the gauss point data to a plane. 

3. The inverse mapping procedure outlined above is then used to map extrapolated 
stresses from the old nodes to the new global gauss point locations saved in step (1). 

4.4 Element Variables Mapped 

In addition to the stresses (<r„, a >y , <j zz ), several other gauss point state variables are 
mapped after remeshing: 
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YieldStress - this value is the current radius of the von Mises yield surface. 

TotalWork - this is a numerical evaluation of the work integral using the trapezoidal 
rule: 


W, = f 'ode (6) 

' Jo 

Plast icWor k - similarly, this is the plastic component of work: 

W r = j I’ede, (7) 

Two of the element gauss point variables are calculated rather than mapped. These are 
logical flags indicating yield state at the gauss point: 

EverYield - indicates whether the gauss point has ever yielded as part of the analysis. 

IsYielding - indicates whether the gauss point is currently yielding. 

Since these two state variables have discrete, logical values, they cannot be mapped. At 
each gauss point they are reevaluated from the current mapped stress state, the current 
yield surface radius, and initial yield surface radius. Using the nomenclature from above: 

EverYield = (YieldStress/YSinit) >=1.02 

IsYielding = (SigE/YieldStress) >0.99 

Where YSinit is the initial radius of the yield surface before hardening, and SigE is 
the von Mises stress evaluated at the current gauss point stress state. 

4.5 Accuracy of the Mapped Data 

There are a number of factors that influence the accuracy of the mapped data. The pri- 
mary factor is the mesh quality. Element quantities are extrapolated to nodes, but ele- 
ment contributions at a node are not averaged. Thus, stress discontinuities across element 
boundaries may exist as a result of inadequate mesh density. Averaged nodal values 
might tend to average out local stress discontinuities; however, local stress gradients may 
be misrepresented. Given the importance of mesh quality on the overall solution, this is 
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probably not a primary concern. That is, if the mesh is not fine enough to accurately 
capture stresses across element boundaries, errors in the mapping results are only con- 
tributing to a much larger problem. 

Another possible error in the mapped data is the relationship between the stress state and 
the yield surface consistency condition for a gauss point that is currently yielding. The 
yielding flags are discrete values that cannot be interpolated, so it is not possible to de- 
termine if a point is yielding without assuming that the mapped stress state and yield sur- 
face radius are consistent. Since the EverYield and IsYielding flags are only used 
for displaying purposes, they are evaluated from the mapped state directly, assuming that 
the mapped state is consistent. 

No special measures are taken to check for the consistency condition or for possible er- 
rors in the stress-strain compatibility other than to re-establish equilibrium in the new 
configuration after the mapping has completed. Typically, this new configuration will 
have new surface area after fracture propagation. Since equilibrium is always re- 
established at the current load and displacement level, this results in stress reduction and 
redistribution as a result of the compliance change in the structure, and a somewhat dif- 
ferent stress state than before the mapping. Numerical experiments have shown that with 
a quality mesh, no problems occur. Habraken and Cescotto (1990) report similar find- 
ings. 

Figure 4.1 through Figure 4.3 provide a qualitative assessment of the mapping errors by 
comparing color contours of stress. The (a) figure shows the stress component and de- 
formed mesh just before crack propagation. The (b) figure shows the stress component 
and the same view of the mesh just after mapping, but before equilibrium has been re- 
stored. The differences in the stress are a direct result of the remeshing and mapping 
process. The figures show that overall quality of the mapped data is good. The largest 
differences appear to be in the vicinity of the original crack tip location (a). 

4.6 Algorithm Efficiency 

Other than solving the nonlinear equilibrium equations, the mapping algorithm is poten- 
tially the most computationally intensive algorithm. For simplicity sake, the initial im- 
plementation mapped the response for all nodes and elements. However, the efficiency 
suffered. For a typical tearing analysis slightly less than V 3 of the total analysis time was 
spent in the mapping algorithm. 
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Consequently, the algorithm was modified to only map for the new nodes and elements 
created during the remeshing process rather than for the entire mesh. Data for other 
nodes and elements is preserved exactly. Further, the algorithm was modified such that 
the search set in the old mesh was limited only to nodes and elements in the remesh re- 
gion. The result is very good performance. For a problem with 6801 nodes and 3288 
elements, the nonlinear solution time for the first crack propagation step is about 1:59 
(mm:ss). Table 1 shows remeshing and mapping times for the original Full mapping al- 
gorithm, the Intermediate algorithm that mapped only for new nodes and elements, and 
the final Optimized algorithm that additionally limits the search set to the nodes and ele- 
ments in the mapping region. 

The number of nodes mapped in this example was 40, and the number of elements 
mapped was 58. It is clear that the algorithmic changes had a dramatic effect on the effi- 
ciency of the algorithm. Since the mapping time is now on the same order as the 
remeshing time and is a small fraction of the total analysis time, the algorithm is deemed 
efficient enough. 

Table 4.1: Comparison of Mapping Algorithm Times (Analysis time: 1:59) 


Algorithm 

Mapping time 

Mapping/Analysis 

Full 

1:51 

93.3% 

Intermediate 

0:09 1 

7.6% 

Optimized 

0:03 

2.5% 


A second order effect on algorithm run-time is the search algorithm for points contained 
in an element. When nodes are mapped from the old mesh to the new mesh, the first step 
is to locate the element in the old mesh that contains a point from the new mesh. A num- 
ber of algorithms are available for solving this problem (Lim, 1992); however, the algo- 
rithm implemented is quite simple. First we search linearly through the list of elements 
and check whether the point lies within the bounding box of each element. For points 
within the bounding box of element, the scan line algorithm is used to determine if the 
point lies within the element itself. 

One additional improvement that was not implemented may improve the mapping effi- 
ciency further. When mapping element gauss point data, each element in the new mesh 
contains several gauss points in close proximity. Rather than searching the entire old 
mesh for each of the gauss points in element from the new mesh, the search could be per- 
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formed only for the first gauss point, then for subsequent gauss points local topology of 
the old mesh can be used for local search, rather than searching the entire mesh. This 
would exploit the fact that all of the gauss points for an element in the new mesh will 
probably lie within the boundaries of a single element in the old mesh, or possibly within 
a small group of elements in the old mesh. 

4.7 Summary 

State variable mapping provides a means for maintaining the plastic history of an elastic- 
plastic analysis as the mesh changes during crack propagation. The inverse isoparametric 
mapping algorithm has proven robust, accurate, and efficient. 
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Figure 4 1 ' Qualitative mapping error assessment. Component: Maximum principal stress, (a) 
just before crack propagation, (b) just after crack propagation, remeshing and mapping, but be- 
fore equilibrium has been restored. 



Figure 4.2- Qualitative mapping error assessment. Component: Effective (von Mises) stress, (a) 
just before crack propagation, (b) just after crack propagation, remeshing and mapping, but be- 
fore equilibrium has been restored. 
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Figure 4.3: Qualitative mapping error assessment. Component: Shear stress, (a) just before 
crack propagation, (b) just after crack propagation, remeshing and mapping, but before equilib- 
rium has been restored. 
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Chapters 


Fracture Analyses Using Nodal Release 


Two fundamentally different algorithms are available for elastic-plastic fracture simula- 
tions in FRANC2D/L: one is designed for efficiency when the fracture occurs along a 
plane of symmetry in the model, and the other is designed for generality when fracture is 
not on planes of symmetry and may involve curving cracks. This chapter presents analy- 
sis results for the former of these two and compares these results with experimental data 
and three-dimensional analysis results. 

The first section of this chapter describes some background information for the analysis 
such as material data, a simple mesh convergence study, and information about how 
three-dimensional analysis is used as a baseline for comparison. The remainder of this 
chapter is devoted to residual strength analyses and their comparison with experimental 
data and three-dimensional analyses. 

Two standard fracture specimens are considered: the M(T) (middle crack tension) and the 
C(T) (compact tension) specimens. These specimens were considered primarily because 
of availability of test data and past analysis experience. First, the plane strain core (PSC) 
height necessary for an accurate residual strength simulation of a 6.0 inch wide C(T) 
specimen is determined. Second, the PSC is verified by comparing analysis predictions 
with experimental data for a 24 inch wide M(T) specimen. This value of the plane strain 
core was then used in all subsequent analyses. Third, the software is used for analysis 
predictions of residual strength for both C(T) and M(T) specimens to show a degree of 
specimen independence, first for configurations with a constant crack-length to speci- 
men-width, then for configurations with constant specimen width and different crack 
lengths. Next, a series of Multi-Site Damage (MSD) calculations demonstrate multiple 
growing cracks with linkup, and show the effects of MSD on residual strength. Finally, 
the effects of specimen thickness on core height is explored. 
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5.1 Analysis Background 

The fracture analyses presented here use the CTOA fracture criterion. All of the simula- 
tions are Mode-I only; a direction criterion is not necessary because symmetry is ex- 
ploited and the nodal release algorithm is used. Figure 5.1 is a schematic of the M(T) test 
specimen. Only % of the specimen need be simulated (symmetry along the centerline and 
fracture surface). The fracture surface symmetry is maintained by the nodal release algo- 
rithm. Figure 5.2 is a schematic of the C(T) test specimen. Similar to the M(T) speci- 
men, the nodal release algorithm maintains the fracture surface symmetry boundary con- 
dition. 

5.1.1 Analysis Material Parameters 

For all of the analyses in this chapter the material considered is 2024-T3 with cracks ori- 
ented in the LT direction (loading parallel to the rolling direction, cracking perpendicular 
to the rolling direction). Table 1 shows the analysis material properties for all analyses in 
this chapter (unless otherwise specified), and Table 2 shows the uniaxial stress-strain data 
for 2024-T3 used in the von Mises yielding algorithm. The critical CTOA value was de- 
termined by Dawicke and Newman (1997) by using three-dimensional finite element 
analysis and comparing failure analysis results with test results. This critical value is 
within the scatter band of experimental CTOA data for 2024-T3. 


Table 5.1: Analysis Material Properties 


Young’s Modulus 

10350 Ksi 

Poisson’s Ratio 

0.3 

Thickness 

0.09 inch 

Yield Stress 

50 Ksi 1 

CTO Ac @ 0.04" 

5.25 deg 


5.1 .2 Mesh Size Convergence Study for Tearing 

The CTOA algorithm implemented in FRANC2D/L allows the crack to advance one or 
more element lengths when the angle made by points on the upper and lower crack sur- 
faces at a distance d = 0.04 inch behind the crack tip reached a critical value. The 0.04 
inch distance was selected based on previous analysis experience (Newman et al., 1993b) 
and on mesh convergence studies (Dawicke, 1996). Figure 5.3 shows results of a simple 
FRANC2D/L convergence study, where the number of elements along the crack face was 
increased, while the critical CTOA was measured at a fixed distance behind the crack tip 
(d = 0.04 inch). The results are compared with similar results from a ZIP2D convergence 
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study (Dawicke, 1996) (the FRANC2D/L results are for linear strain six-noded iso- 
parametric elements (T6); the ZIP2D results are for constant strain triangles (CST)). The 
variation in FRANC2D/L results is less than 1% for all three-mesh densities, thus a 
minimum element size of 0.04 inch on the crack face is deemed acceptable. 


Table 5.2: Uniaxial stress-strain yielding curve 


Stress (Ksi) 

Strain 

0.0 

0.0 

50.0 

0.00483 

56.5 

0.015 

62.5 

0.04 

68.5 

0.1 

71.0 

0.16 

71.0 

0.2 


5.1.3 Three-Dimensional Analyses 

While a wealth of 2024-T3 data is available for M(T) and C(T) specimens, most of this 
data is for only one thickness, and only for limited initial a/W crack ratios. Dawicke and 
Newman (1998) have shown that three-dimensional analysis is capable of accurately rep- 
resenting fracture behavior of thin-sheet material for a wide range of specimen sizes and 
for more than one specimen type. So in addition to experimental data, three-dimensional 
analyses are used as a baseline for comparison with the two-dimensional plane stress, 
plane strain, and plane strain core analyses produced with FRANC2D/L. The three- 
dimensional code is ZIP3D (Shivakumar and Newman, 1990), and Dr. David Dawicke at 
NASA Langley Research Center provided the three-dimensional results (Dawicke, 1997). 

5.2 Constraint Model 

Two advantages of two-dimensional analysis over full three-dimensional analysis are ef- 
ficiency and simplicity. However, with simplicity comes the inability to capture certain 
behavior such as variable through-thickness stress triaxiality, which can be modeled us- 
ing three-dimensional analyses. Stress triaxiality is a significant issue at the tip of a crack, 
even for thin sheet material. This stress triaxiality, or constraint, has received much at- 
tention, but still few models exist to describe constraint in a two-dimensional analysis. 
As previously discussed, plane stress has no constraint, and plane strain introduces too 
much constraint, allowing the plane strain triaxiality to extend too far away from the 
crack tip. Newman et al. (1988) modeled constraint as a simple mixed state of stress with 
plane strain elements near the tip and plane stress elements away from the tip. Figure 5.4 
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is a schematic of the plane strain core (PSC) concept illustrating Newman’s mixed state 
of stress. Analyses in this chapter denoted PSC core use the mixed state of stress and in- 
clude an indication of the core height. 

5.2.1 Plane Strain Core Height 

Different material thicknesses and different test specimens have different levels of con- 
straint, and the height of the PSC can be adjusted to introduce an appropriate level of 
constraint for a given material, thickness and test specimen. For the purposes of this 
chapter, the PSC height was selected by matching the failure load of the two-dimensional 
PSC analysis with test results for 6.0 inch wide C(T) specimens (crack length of a/W = 
0.4). Analyses were run for PSC heights of 0.02, 0.04 and 0.06 inches. Using linear in- 
terpolation, the core height that matched the test results was 0.042 inches. Figure 5.5 
shows a comparison of plane stress, plane strain, and plane strain core results with the 
experimental and three-dimensional results for load against crack extension. The plane 
stress results over-estimated and plane strain results under-estimated the failure load. In 
addition, the crack extension at failure is over-estimated by plane stress and under- 
estimated by plane strain. The plane strain core results closely match the three- 
dimensional results, and represent a good approximation to the test results. 

Similar analyses were run for 12.0 inch wide M(T) specimens with a crack length of 
2a/W =1/3 and plane strain core heights of 0.02 and 0.04 inches. Also using linear inter- 
polation, the core height that matched the three-dimensional M(T) results was 0.034 
inches. 

5.3 Geometry and Specimen Independence 

A necessary condition for a valid fracture criterion is that it be applicable independent of 
specimen type. The critical CTOA fracture criterion is relatively independent of test 
specimen configuration and crack configuration for constant thickness test specimens 
(Dawicke and Newman, 1998). Since the PSC is an analysis mechanism introduced to 
approximate crack-tip constraint, it too must be independent of specimen type and con- 
figuration to be widely applied as a valid analysis approximation. In this section, the 
PSC height determined for the 0.09” thick, 6” wide C(T) specimen is used to predict the 
behavior of other C(T) and M(T) configurations. All analyses were run using a plane 
strain core height of 0.04 inches for convenience. A single core height demonstrates that 
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a core height tuned with one laboratory specimen, can be used to accurately predict the 
behavior of other specimen configuration and other specimen types. 

Figure 5.6 shows a comparison of plane stress, plane strain, and plane strain core results 
with the experimental and three-dimensional results for load against crack extension for a 
24.0 inch wide M(T) specimen. As with the 6.0 inch wide C(T) specimen, the plane 
stress results over-estimated and plane strain results under-estimated the failure load. 
The plane strain core results more accurately model the fracture behavior, over predicting 
the experimental average by about 7%. 

A valid fracture criterion must be independent of specimen configuration. Figure 5.7 
shows predictions of failure load against specimen width for two-dimensional PSC analy- 
ses compared to three-dimensional analysis results for C(T) specimens with constant 
thickness and a/W ratio. The plane strain core analyses accurately predicted the failure 
load for the three largest specimens. A 15% difference in the failure load was observed 
between the two-dimensional and three-dimensional analyses for the 2.0 inch wide C(T) 
specimen. Results indicate the failure of the two inch C(T) specimen is controlled by net 
section yielding, rather than controlled fracture. 

Figure 5.8 shows failure stress against specimen width for M(T) specimens with constant 
thickness and a/W ratio. The results show that the plane strain core can predict failure 
with reasonable accuracy for these specimens as well. For the smallest specimen width, 
the PSC results under-predict the three-dimensional analysis results by about 3%, and for 
the largest specimen width they over-predict by about 3%. 

Figure 5.9 shows predicted failure load against initial crack length for C(T) specimens 
with constant thickness and width. The PSC results agree with the three-dimensional cal- 
culations for all crack lengths. Of the three a/W ratios considered, the PSC results for 
a/W = 0.3 over-estimate the three-dimensional results by 3%; however, for the other two 
ratios the PSC results are nearly identical with the three-dimensional results. Figure 5.10 
shows the failure stress against initial crack length for 24.0 inch wide M(T) specimens 
with constant thickness and width. Again, the results are nearly identical with the three- 
dimensional results. 

These results indicate that the PSC concept can be used to accurately model the fracture 
behavior of both C(T) and M(T) specimens over a range of specimen sizes and initial 
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crack lengths. The largest errors in predicted failure load occurred for the 2.0 inch wide 
C(T) specimen. The FRANC2D/L PSC predictions for the M(T) specimen were within 
2.5% of the ZIP3D predictions of failure stress. As shown in Figure 5.8, the PSC predic- 
tions were lower than the three-dimensional predictions for the smaller specimen widths. 
Conversely, for larger specimen widths the PSC predictions were higher than the three- 
dimensional predictions. 

5.4 Multiple-Site Damage 

Multiple-site damage (MSD) cracking involves two or more cracks that can ultimately 
interact to influence each other's behavior and the behavior of the overall structure. MSD 
can significantly reduce the residual strength of a structure. The ability of software to 
predict residual strength in the presence of MSD is particularly important. Figure 5.1 1 is 
a schematic of the idealized geometry for eight different cracking configurations 
(Dawicke and Newman, 1998). The specimen was a standard 12.0 inch wide M(T) con- 
figuration, but with different cracking configurations. The actual crack configurations of 
the specimens were not exactly symmetric, so the ligament lengths and crack lengths 
were averaged between the two sides to obtain these symmetric configurations for the 
analysis. The first three configurations have single center cracks with 2a/W ratios of 
4/12, 5/12 and 6/12, respectively. The subsequent configurations are multiple crack con- 
figurations. The material was 0.09 inch thick 2024-T3, and all cracks were fatigue pre- 
cracked at low stress levels. 

The same material data used in previous analyses was used for each of the analyses of the 
eight configurations. As before the plane strain core height was 0.04 inch and the critical 
CTO A was 5.25°. The analyses were run under displacement control with symmetry 
boundary conditions for both the specimen load center-line and for the crack face. 

Table 3 summarizes the comparison of experimental link-up and failure loads with the 
predicted FRANC2D/L results. FRANC2D/L predicts the failure load within 6%, and 
predicts link-up load within 2.4% for all of the cases. Configurations G and H failed at 
link-up. For these cases, FRANC2D/L correctly predicted link-up before failure. 


Table 5.3: Comparison of FRANC2D/L failure and link-up stresses with test results. 


Cfg 

Exp 

Failure 

(ksi) 

FRANC2D/L 

Failure 

(ksi) 

Failure 

Error 

Exp 

Link-Up 

(ksi) 

FRANC2D/L 

Link-Up 

(ksi) 

Link-Up 

Error 

A 

31.3 

30.6 
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B 

27.7 

26.7 

-2.4% 




C 

23.3 

22.4 

-4.0% 




D 

30.3 

30.6 

+1.0% 

13.2 

13.1 

-1.0% 

E 

28.2 

28.1 

-0.4% 

21.6 

21.1 

-2.4% 

F 

23.0 

22.4 

-2.7% 

19.7 

19.8 

+0.5% 

G 

37.5 

36.9 

-1.6% 

At Failure 

19.6 


H 

23.4 

24.8 

+6.0% 

At Failure 

22.1 & 17.2 



Several interesting observations can be made through comparisons of the various con- 
figurations. First, by comparing configurations A and D, we can see that the remaining 
ligament of the two-crack configuration has little effect on the residual strength, reducing 
the strength only slightly. Comparing configuration A with configurations F and H, we 
see that in both cases the small damage in front of the main crack significantly reduces 
the strength of the specimen. Comparing configuration C with F we see that F behaves 
more like a single 6” crack, rather than a 4” crack with damage. That is, the damage 
makes the effective length look like the total of the main crack plus damage. Configura- 
tion G, which had the highest strength, shows that damage in front of a main crack will 
not always reduce the strength to the level of the total length including the damage. 
However, assuming the reduction occurs will be conservative. 

Figure 5.12 through Figure 5.14show the crack tip link up process for Configuration F, 
which has a 4 inch lead crack with a single pair of MSD cracks in front of the lead crack. 
The MSD cracks are each about 0.56" long and the remaining ligament between each 
MSD crack and the lead crack is about 0.48" long. The figures all show color contours of 
the von Mises stress using the same stress range and deformed geometry magnified xl5. 
Figure 5.12 shows the condition just before Tip 1 is ready to propagate. The fracture 
criterion at Tips 2 and 3 is slightly below critical. Figure 5.13 shows the condition where 
the remaining ligament between Tips 1 and 2 is ready to break. Tip 2 is at the fracture 
criterion. As the crack face forces are relaxed for Tip 2, load redistribution causes Tip 1 
to exceed the fracture criterion. Crack face forces are subsequently relaxed for both Tip 1 
and Tip 2. The resulting load redistribution causes Tip 3 to meet the fracture criterion, 
and crack face load relaxation ensues at Tip 3. Each node along the crack face has a 
count of the number of steps required for load relaxation. Only after all nodes along the 
crack face are finish relaxing does the far field boundary condition start increment again. 
Essentially, all cracks must be stable before the boundary condition increments will re- 
sume. Figure 5.14 shows the next sable condition after Tips 1 and 2 have merge, break- 
ing the ligament between the lead crack and the MSD crack. The model now contains a 
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single crack. Figure 5.14 shows clearly residual stresses in the wake of the ligament, and 
the effect of the residual stresses on the deformed shape. 

5.5 Effects of Material Thickness on Plane Strain Core Height 

The fact that material thickness affects constraint has been understood for some time. 
ASTM standards specify requirements that must be met for plane strain fracture tough- 
ness to determine a Ki c critical stress-intensity factor. These standards are, in effect, es- 
tablishing that enough constraint must exist to consider the test sufficiently close to plane 
strain conditions. Thickness will have a similar effect for elastic-plastic fracture. This 
section presents elastic-plastic results used to explore how plane strain core height varies 
with the thickness of a specimen. 

Just as the critical stress-intensity factor varies with thickness, the critical CTOA varies 
with thickness. However, no significant studies exist to characterize CTOA as a function 
of thickness for 2024-T3. As a result, for the purposes of this study, CTOA is assumed 
constant at 5.25° for all thicknesses analyzed. 

Figure 5.15 and Figure 5.16 show the variation in failure load as a function of PSC height 
for a 6” C(T) specimen and 1 2” M(T) specimen, respectively. Three-dimensional analy- 
sis results were used to find the “correct” core height to predict failure for the respective 
specimen and thickness. This same process was repeated for both the 6” C(T) and 12” 
M(T) specimens over a range of thicknesses from 0.025” to 0.25”. The experimentally 
verified critical CTOA for the B=0.09” thick specimen was used for the full range of 
specimens for the study. While in reality the critical CTOA will almost certainly vary 
with thickness for these specimens, the qualitative results are important. Figure 5.17 
shows how the PSC height varies with specimen thickness. The graph shows that the re- 
lationship between core height and thickness is not linear. For all thicknesses below 
about 0.18” the C(T) core height is larger than the M(T) core height. A transition occurs 
at about 0.018” after which the M(T) core height is larger then the C(T) core height. The 
general trend in the curves for both specimens is about the same. 

The implication of the plateau for the smaller thicknesses may be that even for very thin 
specimens, some level of stress traxiality exists at the crack tip, resulting in constraint 
that must be modeled for the most accurate residual strength analyses. 
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5.6 Summary and Conclusions 

This chapter presented analysis results for plane stress, plane strain, and plane strain core 
approximations. Three-dimensional analyses were used in conjunction with experimental 
data as a basis for comparison. Plane stress and plane strain results do not capture the 
variable constraint that exists throughout the specimen during fracture. The plane strain 
core was able to characterize constraint and produced more accurate residual strength 
analyses. The plane strain core proved effective for two different standard laboratory 
specimen types (M(T) and C(T)), for a variety of specimen sizes, and for a variety of ini- 
tial crack-length to width ratios as well as MSD cracking configurations. 

The conclusions are: 

• The plasticity implementation in FRANC2D/L is accurate and effective for Mises 
type materials. 

• The plane strain core concept can approximate the three-dimensional constraint ef- 
fects present during fracture of thin-sheet aluminum. 

• A single value of CTOA determined from tests on a C(T) laboratory specimen were 
used to predict failure for other C(T) geometries and for M(T) specimen geometries, 
including variations in both specimen size and initial crack length. 

• The plane strain core concept can be used to accurately simulate fracture and link-up 
for MSD cracking configurations. 

• The relationship between material thickness and plane strain core height is not linear. 

• Even the thinnest specimen considered (0.025”) required a PSC to best model the ref- 
erence three-dimensional results. 
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Figure 5.1: Schematic of a deformed middle-crack tension specimen showing how symmetry is 
exploited. 
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Figure 5.2: Schematic of a deformed compact tension specimen showing how symmetry is ex- 
ploited. 



Figure 5.3: Mesh convergence study results for the 3" M(T) specimen. The small strain, six- 
noded triangle implemented in FRANC2D/L performs well at each of the element discretizations. 
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Figure 5.4: Plane strain core concept. 



Figure 5.5: Comparison of simulated load vs. crack extension with experimental data for 6.0 inch 
wide C(T) specimens. The plane strain core closely matches the 3D result, and provides a good 
approximation to the experimental results. 
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Figure 5.6: Comparison of simulated load vs. crack extension with experimental data for 24.0 
inch wide M(T) specimens. The plane strain core provides a much better approximation to the 
experimental results than either plane stress or plane strain. 



Width (In) 


Figure 5.7: Predicted failure load vs. specimen width for C(T) specimens with constant thickness 
and a/W. 
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Figure 5.8: Predicted failure stress vs. specimen width for M(T) specimens with constant thick- 
ness and 2aA N. 



Figure 5.9: Predicted failure load vs. initial crack length for C(T) with constant 
thickness and width. 
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Figure 5.10: Predicted failure stress vs. initial crack length for M(T) with constant thickness and 
width. 
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Figure 5.11: Additional single crack and MSD configurations analyzed (all configurations are 
symmetric). 
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Figure 5.12: Impending crack growth at Tip 1 for MSD Configuration F. Color contours of von 
Mises stress (Ksi) on deformed body magnified x15. 
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Figure 5.13: Remaining ligament between the lead crack and MSD crack ready to break. Color 
contours of von Mises stress (Ksi) on deformed body magnified x15. 
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Figure 5.14: Lead crack and MSD crack have merged, leaving a single crack. Color contours of 
von Mises stress (Ksi) on deformed body magnified x15. 
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Figure 5.1 5: Variation in failure load as a function of PSC height for a 6” C(T) specimen . 
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Figure 5.16: Variation in failure load as a function of PSC height for a 12" M(T) specimen. 
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Figure 5.17: Variation of plane strain core height with specimen thickness. 
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Chapter 6 


Fracture Analyses with Remeshing 


Modeling curvilinear crack growth without knowing the crack path in advance requires 
mechanisms to predict and represent an evolving crack geometry. In this chapter we pre- 
sent a series of analyses to demonstrate the remeshing approach to elastic-plastic crack 
growth. First, analyses that use the remeshing approach will be compared to analyses 
using the nodal release approach. These analyses are Mode-I only. Next, comparisons 
are made with experimental data for curving cracks. 

6.1 Analysis Background 

Some of the fracture analyses presented here are Mode-I and some are Mode-I/II. How- 
ever, all of the analyses were generated using the remeshing and mapping algorithms and 
the direction and fracture criteria discussed in previous chapters. The fracture criterion 
used is the magnitude of the CTOD vector, and the direction criterion used is the maxi- 
mum tensile stress evaluated at the crack tip node. 

Three different specimens were used for the analyses in this chapter: 6” C(T), 12” M(T) 
and a mixed-mode I/II Arcan specimen. The C(T) and M(T) specimens are Mode-I only 
test specimens. Numerical results from Chapter 5 for the specimens are used for com- 
parison with the results from the remeshing-mapping approach. Tables 5.1 and 5.2 give 
material properties for 2024-T3 aluminum used for the Mode-I comparisons. 

The mixed-mode analyses are simulations of experiments conducted by Amstutz et al. 
(1995a, 1995b). Figure 6.1 is a schematic of the Arcan specimen used by Amstutz for the 
mixed-mode testing. The specimen allows a full range of loading from Mode-I to Mode- 
II. The Amstutz et al. (1995a, 1995b) results are also for 2024-T3 aluminum, so Tables 

5.1 and 5.2 described the material properties for the mixed-mode analyses (same thick- 
ness as for the C(T) and M(T) tests). As for the M(T) and C(T) specimens, all Arcan 
specimens were fatigue pre-cracked at low stress levels in the LT orientation, and all tests 
were conducted under displacement control. 
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The value of critical CTOD used for the mixed-mode analyses deserves special mention. 
As discussed previously, CTOD for 2024-T3 varies with the rolling direction of the mate- 
rial by about 25% (Amstutz et al., 1995b). Dawicke (1996) found CTOA at 0.04” behind 
the crack tip to be 5.25 degrees. The equivalent CTOD is 0.00366". This value of CTOD 
falls in the midrange of the experimental data measured by Amstutz et al. (1995b). The 
most accurate simulation would account for the variation in CTOD as a function of crack 
tip orientation with respect to material rolling direction. However, for this analysis only 
the single value of CTOD based on Dawicke’s result is used. 

6.2 Mode-1 Analyses 

The results of this section serve mainly as verification that the remeshing-mapping ap- 
proach gives the same results as the nodal release approach. Two standard test specimen 
configurations are used for comparisons in this section: compact tension specimens and 
middle crack tension specimens. The results from Chapter 5 for FRANC2D/L where 
C(T) and M(T) specimens were analyzed using nodal release are included for compari- 
sons with the results from this chapter. Two specimen configurations are considered: a 
6” wide C(T) specimen with a/W = 0.4, and a 12” wide M(T) specimen with 2a/W =1/3. 

Figure 6.2 shows the mesh for the C(T) specimen, and Figure 6.3 shows the mesh near 
the crack tip. The crack tip element size is 0.02”, and the mesh into which the crack 
propagates is uniform with a characteristic element size of about 0.02”. Since CTOD is 
measured at 0.04" behind the crack tip, there are two elements along the crack face be- 
hind the tip to the measurement point. 

Figure 6.4 shows the load-crack extension curves for the 6” wide C(T) specimen. The 
solid lines show results for automatic remeshing and dashed lines show results for nodal 
release from Chapter 5. Also included in Figure 6.4 are test data from Dawicke and 
Newman (1998). As expected, although the nodal release and automatic remeshing algo- 
rithms are completely different, the results match well. 

Figure 6.5 shows the mesh for the 12” wide M(T) specimen. The crack tip element size 
is 0.02”, and the mesh into which the crack propagates has a characteristic element size 
of about 0.04”. Figure 6.6 shows the load-crack extension curves for the 12”wide M(T) 
specimen. The solid lines show results for automatic remeshing and dashed lines show 
results for nodal release from Chapter 5. Also included in Figure 6.6 are test data from 
Dawicke and Newman (1998). The results from the current analysis using automatic 
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propagation with remeshing are indicated by the dashed lines from a plane stress case. 
Again, the results match well. 

6.3 Mixed-Mode I/ll Analyses 

Figure 6.1 is a schematic showing the definition of loading angle and specimen dimen- 
sions for the Arcan specimen and fixture tested by Amstutz et al. (1995a, 1995b). Single 
point bolt loads can be applied on opposite sides of the fixture to create loading that is 
pure Mode-I, pure Mode-II, or one of five intermediate values. Figure 6.7 shows a mesh 
of the Arcan specimen for a 30° loading angle. To simplify post processing, the mesh 
was rotated for each loading angle and loading always applied vertically. The crack tip 
element size is 0.02” and the surrounding element size is about 0.04”. Loading angles 
from 0° to 90° in 15° increments were analyzed. 

The analyses were performed by loading the specimen vertically under displacement 
control. When the CTOD reached the critical value, the crack propagated in the direction 
of maximum tensile stress. As shown in Figure 6.8, significant plasticity occurs in the 
specimen. Thus, propagation is occurring through material that has experienced signifi- 
cant yielding. 

Figure 6.9 compares the crack path recorded during experiments and the crack path pre- 
dicted during the FRANC2D/L analysis for loading angles of 0, 30, and 60 degrees. The 
maximum tensile stress and CTOD fracture criterion predict the crack path well. The 
FRANC2D/L analyses tend to under-predict the initial kinking angle and then slightly 
over-predict the later propagation angle. For loading angles greater than 60° the material 
fails in shear. As discussed in Chapter 3, the maximum tensile stress criterion cannot 
predict shear failure, so these analyses are not included. 

The analyses also predict the residual strength of the specimen for the different loading 
angles. Figure 6.10 compares the experimental load-crack extension with the 

FRANC2D/L results for 0 and 60 degree loading angles. In general terms, the 
FRANC2D/L analyses predict the trend of lower strength for the 60 degree case, com- 
pared to the 0 degree case. In addition, the analyses predict the cross-over in load-crack 
extension curves displayed by the experimental data. More specifically, the 
FRANC2D/L result over predicts the experimental results by 9% for the zero degree case 
and over predicts by 6% for the 60 degree case. 


72 



Figure 6.11 compares the experimental load-crack extension with the FRANC2D/L re- 
sults using plane stress and plane strain assumptions for the 30 degree loading angle. The 
plane stress analysis over-predicts by 10%, while the plane strain analysis significantly 
over-predicts the experimental result. The character of the plane strain result is similar to 
that for the previous M(T) analysis, in that the plane strain curve does not cross the plane 
stress curve. For the previous C(T) result, the curves do cross. 

Results are not presented for the 75 and 90 degree analyses. As previously discussed, 
these loading cases fail in Mode-II dominated propagation. The maximum tensile stress 
criterion is not valid for those cases. 

6.4 Summary 

The results presented show that use of an elastic-plastic material model with a CTOD 
fracture criterion and maximum tensile stress direction criterion can predict the residual 
strength and crack path for test specimens with significantly different fracture behavior. 
The Arcan test results show significant crack curving. The analyses capture this behav- 
ior, even in the presence of significant plasticity. The analyses also predict residual 
strength within an error of approximately 10%. Because the value of critical CTOD used 
in these analyses was obtained independently and not adjusted in these analyses, the re- 
sults demonstrate predictive capability for elastic-plastic fracture. 

Reliable analysis relies on accurate plasticity states for the material into which the crack 
is propagating as well as capturing the residual stresses along the crack face. For the 
model presented here, stresses and the nodal displacements (and hence, strains) are 
mapped from one mesh to another during propagation. Ensuring a quality mesh and en- 
forcing equilibrium after each propagation step was adequate to obtain reasonable frac- 
ture predictions. Similarly, the results show that reasonable fracture predictions can be 
obtained using a well refined mesh, without resorting to special elements at the crack tip. 

It is clear from the results that a plane strain core model to capture crack tip constraint 
improves the accuracy of the predicted failure loads for the M(T) and C(T) specimens. It 
is not clear how well this approach will work for a Mixed-Mode growing crack. Equally 
important is the unknown effect of a constraint model on the criterion used for direction 
prediction. These topics are areas for future research. 
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Figure 6. 1 : Schematic of Arcan test specimen and fixture. 



Figure 6.2: Finite element mesh for C(T) analysis. 
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Figure 6.3: Detail of mesh near crack tip for 6” wide C(T) specimen. 



Figure 6.4: Load crack extension curve for C(T) analysis. 
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Figure 6.5: Finite element mesh for 12" wide M(T) analysis. 
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Figure 6.7: Finite element mesh for Arcan analyses (30 degree load angle). 
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Figure 6.8: Shaded region showing elements that have yielded after 10 mm of crack propagation 
(30 degree load angle). 
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Figure 6.9: Observed and predicted crack paths for Arcan specimen. 



Figure 6.10: Load crack extension curve for Arcan analyses. 
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Figure 6.1 1 : Comparison of plane stress and plane strain analyses for Arcan specimen with 30 
degree load angle. 
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Chapter 7 


Summary 


This chapter presents a summary of this work. Each chapter is summarized, then conclu- 
sions are drawn based on experience with the implementation and on results from chap- 
ters of example applications. Finally, brief recommendations outline areas for near-term 
improvement and long-term research. 

7.1 Chapter Summaries 

Chapter 1 introduced the motivation and context for this work and introduced the crack 
growth methodology, which served as the outline for the remainder of this work. 

In Chapter 2, we discussed fracture and direction criteria that may be appropriate for 
Mixed-Mode fracture. Energy based fracture integrals have been used extensively in the 
past in the form of resistance curves, and will continue to be used in the future for known 
configurations, but their applicability for Mixed-Mode fracture, especially in multiple 
cracks scenarios, is in question. For this work the fracture criterion implemented was the 
critical crack tip opening displacement. While this criterion also has limitations and has 
not been proven widely, it does perform extraordinarily well in the circumstances tested. 
Similarly, no single theory exists to predict the crack direction. Recent results from the 
literature indicate that a combined theory of maximum tensile stress and maximum shear 
stress may predict direction for materials that fail in tension under some loading condi- 
tions and in shear under other loading conditions. For this work the maximum tensile 
stress criterion was implemented. While this criterion will only predict the crack growth 
direction correctly for materials with loading that causes tensile dominated fracture, im- 
plementation was relatively straightforward and has proven robust. Conveniently, for 
this criterion no special material dependent properties are needed. 

Chapter 3 discussed various approaches to elastic-plastic fracture from the literature and 
describes the two different approaches to fracture implemented: nodal release and auto- 
matic remeshing with mapping. The algorithms have been implemented to model elastic- 
plastic fracture under a number of different circumstances. The nodal release algorithm 
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is specialized to work for either of two different situations: cracking along an edge of 
symmetry or cracking along a predefined curving crack path. The automatic remeshing 
algorithm can be used either in a step-by-step manual mode or in fully automatic mode. 
In fully automatic mode the direction is predicted from maximum tensile stress at the 
crack tip. Crack tip constraint can be modeled during nodal release using the plane strain 
core. Specialized algorithms allow multi-site damage during nodal release. Manual 
remeshing and mapping can be used to merge MSD crack tips when automatic propaga- 
tion is interrupted by crack tips in close proximity to one another. 

In Chapter 4, we discuss the state variable mapping algorithm that maintains the element 
and nodal solution when remeshing is used to propagate cracks. The inverse iso- 
parametric mapping algorithm has proven robust, accurate, and efficient. The algorithm 
uses knowledge of the remeshing zone around each crack tip and only performs mapping 
for these nodes and elements. Test problems show that the mapping time is approxi- 
mately the same as the remesh time, which is a small fraction of the total analysis time. 

Chapter 5 presented analysis results for two different symmetry specimen types: compact 
tension and middle crack tension. Plane stress and plane strain results do not capture the 
variable constraint that exists throughout the specimen during fracture. The plane strain 
core implementation was able to characterize constraint and produced more accurate re- 
sidual strength analyses. A simple study on constant thickness configurations showed 
that a single plane strain core height determined by strength comparisons with experi- 
mental and three-dimensional data can be used to predict the behavior of a wide variety 
of other specimen types and configurations. In addition, multi-site damage simulations 
were accurate for both residual strength and link-up loads when the PSC was included. A 
second study considered a single configuration for the M(T) and C(T) specimens and 
varied the specimen thickness. For each thickness of each specimen, the best plane strain 
core height was found by matching failure load with three-dimensional finite element re- 
sults. The plane strain core height was not the linear function of thickness, and a plane 
strain core was necessary even for the thinnest specimen, indicating that crack tip con- 
straint is important even for very thin specimens. 

Chapter 6 presented Mixed-Mode analysis results. The results presented show that use of 
an elastic-plastic material model with a CTOD fracture criterion and maximum tensile 
stress direction criterion can predict the residual strength and crack path for test speci- 
mens with significantly different fracture behavior. The Arcan test results show signifi- 
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cant crack curving. The analyses capture this behavior, even in the presence of signifi- 
cant plasticity. The analyses also predict residual strength within an error of approxi- 
mately 10%. Similarly, the results show that reasonable fracture predictions can be ob- 
tained using a well refined mesh, without resorting to special elements at the crack tip. 
Because the value of critical CTOD used in these analyses was obtained independently 
and not adjusted in these analyses, the results demonstrate predictive capability for elas- 
tic-plastic fracture. It is clear from the results that a plane strain core model to capture 
crack tip constraint improves the accuracy of the predicted failure loads for the M(T) and 
C(T) specimens. It is not clear how well this approach will work for a Mixed-Mode 
growing crack. Equally important is the unknown effect of a constraint model on the 
criterion used for direction prediction. 

7.2 Conclusions 

The primary result of this work is an elastic-plastic fracture simulation tool implemented 
in a relatively general way, such that as new technologies emerge (e.g. fracture or correc- 
tion criteria), they can be implemented with relative ease. In addition, the software is im- 
plemented within a relatively user-friendly framework, such that researchers and prac- 
ticing engineers can apply the model to realistic problems. 

The comparisons of analysis results with experimental data show that the model can pro- 
duce accurate residual strength and crack path predictions for a variety of specimens and 
crack configurations. The critical crack tip opening displacement vector proved effective 
as a fracture criterion for Mixed-Mode analyses. While the maximum tensile stress di- 
rection criterion is clearly not applicable in all situations, its simplicity makes it attractive 
as a general purpose direction criterion for materials that fail under tensile dominated 
loading conditions. The results also show that a two-dimensional approximation needs 
some form of constraint model to capture the stress triaxiality near the crack tip and pro- 
duce the most accurate strength predictions. 

7.3 Recommendations for Future Work 

In the near-term several improvements to the software would prove useful. First, a plane 
strain core implementation to capture some measure of constraint during remeshing- 
mapping analyses could increase the accuracy of residual strength predictions. Care 
should be taken to investigate the effect of having elements that change from plane stress 
to plane strain as the core extends. Also, the effects of the core on the local tip stress 
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state should be investigated to ensure that no artificial influence on the direction is intro- 
duced. Second, recent advances in criteria that can predict crack growth direction under 
both tensile and shear dominated situations should be investigated. Finally, the J and T 
integrals may prove to be useful additions to the software. The J R curve has been a stan- 
dard in the industry for a number of years and will continue to be used in the future under 
many situations. The T R curve may in the future prove to have more utility than the J R 
curve, once its limitations are better understood. 

Another improvement to the software would have long-term benefit to the software as 
research tool. A crack-face contact model would open up the possibility for applying 
FRANC2D/L in a variety of circumstances, including: 

• explicit modeling of fatigue crack growth on a cycle by cycle basis, where crack face 
closure and plasticity in the wake is modeled explicitly; 

• softening behavior at the crack tip opens of the possibility of modeling elastic-plastic 
crack growth for materials like concrete, which may have a large process zone; 

• crack formation, where the softening behavior associative with void growth might be 
modeled along the interface. 

A final recommendation for future work is to implement a finite strain model. Many high 
toughness materials exhibit significant blunting with corresponding large strains around a 
crack tip. Computer models of these materials will require a finite strain model to accu- 
rately model fracture. 


84 



References 


Amstutz, B.A., Sutton, M.A., Dawicke, D.S., and Newman, J.C., Jr., 1995a, “An Experi- 
mental Study of CTOD for Mode I/Mode II Stable Crack Growth in Thin 2024-T3 
Aluminum Specimens,” Fracture Mechanics: 26 th Volume, ASTM STP 1256. 

Amstutz, B.A., Sutton, M.A., Dawicke, D.S., and Boone, M.L., 1995b, “Effects of Mixed 
Mode I/II Loading and Grain Orientation on Crack Initiation and Stable Tearing in 
2024-T3 Aluminum,” Fatigue and Fracture Mechanics: 27 th Volume, ASTM STP 
1296. 

Anderson, H., 1973, “A Finite Element Representation of Stable Crack Growth,” Journal 
of the Mechanics and Physics of Solids, Vol 21, pp. 337-356. 

Anderson, T.L., 1995, Fracture Mechanics, Fundamentals and Applications, 2 nd Ed., 
CRC Press, Boca Raton. 

Bakuckas, J.G., Jr., Tan, T.M., Lau, A.C.W., and Awerbuch, J., 1993, “A Numerical 
Model for Predicting Crack Path and Modes of Damage in Unidirectional Metal Ma- 
trix Composites,” Journal of Reinforced Plastics and Composites, Vol. 12, March 
1993. 

Bathe K.J., 1996, Finite Element Procedures, Prentice Hall, Englewood Cliffs, NJ. 

Brown, S.B., and Song, H., 1993, “Rezoning and Dynamic Substructuring Techniques in 
FEM Simulations of Welding Processes,” Journal of Engineering for Industry, Vol. 
115, pp. 415-423, November 1993. 

Camacho, G.T. and Oritz, M., 1996, “Computational Modeling of Impact Damage in 
Brittle Materials,” International Journal of Solids and Structures, Vol. 33, No.20-22, 
pp.2899-2938.. 

Chao, Y.J. and Liu, S., 1998, “On the Failure of Cracks Under Mixed-Mode Loads,” In- 
ternational Journal of Fracture, Vol. 87, No. 3, pp. 201-223. 

Chavez, P.F., 1983, “Automatic Procedures in Evolutionary Finite Element Calculations: 
Restoration of Deteriorated Meshes, Data Transfer Between Meshes, and Mesh Re- 
finement,” Ph.D. Thesis, Cornell University, Ithaca, New York, January 1983. 

Crisfield M.A., 1991, Non-linear Finite Element Analysis of Solids and Structures, 
Wiley, West Sussex, England. 

Dalle Donne, C. and Doker, H., 1997, “Plane Stress Crack Resistance Curves of and In- 
clined Crack under Bi-Axial Loading,” Multi-Axial Fatigue and Deformation Testing 
Techniques, ASTM STP 1280, pp. 243-263. 



Dawicke, D.S. and Newman, J.C., Jr., 1997, “Evaluation of Fracture Parameters for Pre- 
dicting Residual Strength of Multi-site Damage,” First Joint DoD/FAA/NASA con- 
ference on Aging Aircraft, Ogden, Utah, July 8-10. 

Dawicke, D.S. and Newman, J.C., Jr., 1998, “Residual Strength Predictions for Multiple 
Site Damage Cracking Using a Three-Dimensional Finite Element Analysis and a 
CTOA Criterion,” Submitted for Publication in Fatigue and Fracture Mechanics: 

29th Volume, ASTMSTP 1321. 

Dawicke, D.S. and Sutton, M.A., 1993, “Crack Tip Opening Angle Measurements and 
Crack Tunneling Under Stable Tearing in Thin Sheet 2024-T3 Aluminum Alloy,” 
NASA CR- 191 523. 

Dawicke, D.S., 1996, “Residual Strength Predictions Using a CTOA Criterion,” FAA- 
NASA Symposium on Continued Airworthiness of Aircraft Structures, Atlanta, GA, 
August 1996. 

Dawicke, D.S., 1997, Personal communication, NASA Langley Research Center, 
Hampton, Virginia. 

de Koning, A.U., 1977, “A Contribution to the Analysis of Quasi Static Crack Growth in 
the Sheet Materials,” in Fracture 1977, Proceedings of the 4 th International Confer- 
ence on Fracture, Vol. 3, pp. 25-31. 

Demofonti, G. and Rizzi, L., 1991, “Experimental Evaluation of CTOA in Controlling 
Unstable Ductile Fracture Propagation,” ESIS Pub. 9, pp. 693-703. 

Deng, X., 1990, “Dynamic Crack Propagation in Elastic-Plastic Solids,” Ph.D. Thesis, 
California Institute of Technology, Pasadena, CA. 

Erdogan, F. and Sih, G.C., 1963, “On the Crack Extension in Plates Under Plane Loading 
and Transverse Shear,” Journal of Basic Engineering, December, pp. 519-527. 

Espinosa, H.D., Emore, G., and Zavattieri, P., 1996, “Computational Modeling of Geo- 
metric and Material Nonlinearities with an Application to Impact Damage in Brittle 
Materials,” Advances in Failure Mechanisms in Brittle Materials, ASME AMD-Vol. 
219, pp. 119-161. 

Gerstle, W.H., 1986, “Finite and Boundary Element Modeling of Crack Propagation in 
Two- and Three-Dimensions using Interactive Computer Graphics,” Ph.D. Thesis, 
Cornell University, Ithaca, NY. 

Ghosal, A.K. and Narasmhan, R., 1994, “A Finite Element Analysis of mixed-mode 
Fracture Initiation by Ductile Failure Mechanisms,” Journal of Mechanics and Phys- 
ics of Solids, Vol. 42, No. 6, pp. 953-978. 

Gondhalekar, S.R., 1992, “Development of a Software Tool for Crack Propagation analy- 
sis in Two-dimensional Layered Structures,” Master’s Thesis, Department of Me- 
chanical Engineering, Kansas State University, Manhattan, KS. 

Habraken, A.M. and Cescotto, S., 1990, “An Automatic Remeshing Technique for Finite 
Element Simulation of Forming Processes,” International Journal for Numerical 
Methods in Engineering, Vol. 30, pp. 1503-1525. 


86 



Harris, C.E., 1990, “NASA Aircraft Structural Integrity Program,” NASA TM 102637. 

Harris, C.E., 1996, NASA Airframe Structural Integrity Working Group Meeting, NASA 
Langley, Hampton, VA, March 21-22, 1996. 

Hill, R., The Mathematical Theory of Plasticity, Oxford University Press, London, 1950. 

Horn, C.L. and McMeeking, R.M., 1980, “Large Crack Tip Opening in Thin Elastic- 
Plastic Sheets,” International Journal of Fracture, Vol. 45, pp. 103-122. 

Hughes, T.J.R., 1984, “Numerical Implementation of Constitutive Models: Rate- 

independent Deviatoric Plasticity,” Theoretical foundation for large-scale computa- 
tions for nonlinear material behavior (Mechanics of elastic and inelastic solids; 6), 
Nemat-Nasser, Asaro, and Hegemier eds., Martinus Nijhoff Publishers, Dordrecht, 
1984. 

Hutchinson, J.W., 1968, “Singular Behavior at the End of a Tensile Crack Tip in a Hard- 
ening Material,” Journal of the Mechanics and Physics of Solids, 16, pp. 1 3-3 1 , 1 968. 

Hutchinson, J.W. and Paris, P.C., 1979, “Stability analysis of J-Controlled Crack 
Growth,” ASTM STP 668, pp. 37-64. 

Inglis, C.E., 19 13, “Stresses in a Plate Due to the Presence of Cracks and Sharp Comers,” 
Transactions of the Institute of Naval Architects, Vol. 55, pp. 219-241, 1913. 

Irwin, G.R., 1957, “Analysis of Stresses and Strains Near the End of a Crack Traversing a 
Plate,” Journal of Applied Mechanics, Vol. 24, pp. 361-364. 

James, M.A., 1997, “ Residual Strength Calculations for Single and Multiple-Site Dam- 
age Cracks,” First Joint DoD/FAA/NASA conference on Aging Aircraft, Ogden, 

Utah, July 8-10. 

Kato, K., Lee, N.S., and Bathe, K.J., 1993, “Adaptive Finite Element Analysis of Large 
Strain Elastic Response,” Computers & Structures, v 47, n 4/5 pp. 829-855. 

Kfouri, A.P. and Brown, M.W., 1995, “A Fracture Criterion for Cracks under Mixed 
Mode Loading,” Fatigue and Fracture of Engineering Materials and Structures, Vol. 
18, No. 9, pp. 959-969. 

Kfouri, A.P, 1996, “Crack Extension under Mixed-Mode Loading in an Anisotropic 
Mode-Asymmetric Material In Respect of Resistance to Fracture,” Fatigue and 
Fracture of Engineering Materials and Structures, Vol. 19, No. 1, pp. 27-38. 

Krieg R.D. and Key, S.W., 1976, “Implementation of a time independent plasticity theory 
into structural computer programs,” Constitutive Equations in Viscoplasticity: Com- 
putational and Engineering Aspects, AMD-20, ed. J. A. Stricklin et al., ASME, New 
York, pp. 125-138. 

Krieg, R.D. and Krieg D.B., 1977, “Accuracies of numerical solution methods for the 
elastic-perfectly plastic model,” Trans. ASME, pp. 510-515, November. 

Krishnan, S., 1994, “A Finite Element Model of Crack Growth in Layered Structures 
with Bending,” Master’s Thesis, Department of Mechanical Engineering, Kansas 
State University, Manhattan, KS. 


87 



Kuang, J.H. and Chen, Y.C., 1997, “The Tip Plastic Energy Applied to Ductile Fracture 
Initiation under Mixed-Mode Loading,” Engineering Fracture Mechanics, Vol. 58, 
No. 1/2, pp. 61-70. 

Lee, J., Liebowitz, H., Lee, K., 1996, “The Quest of a Universal Fracture Law Governing 
the Process of Slow Crack Growths” Engineering Fracture Mechanics, Vol. 55, No. 

1, pp. 61-83. 

Lee, K., Lee, J., Liebowitz, H., 1997, “Finite Element Analysis of the Slow Crack 

Growth Process in Mixed Mode Fracture,” Engineering Fracture Mechanics, Vol. 56, 
No. 4, pp. 551-557. 

Lim, I.L., Johnston, I.W., Choi, S.K., and Murti, V., 1992, “Improved Numerical Inverse 
Isoparametric Mapping Technique for 2D Mesh Rezoning”, Engineering Fracture 
Mechanics, Vol. 41, No. 3, pp. 417-435, Feb 1992. 

Lim, I.L., 1992, “Fracture Propagation in the Soft Rock,” Ph.D. Thesis, Monash Univer- 
sity, Australia. 

Maccagno, T. and Knott, J., 1989, “The Fracture Behavior of PMMA in Mixed Modes I 
and II,” Engineering Fracture Mechanics, Vol. 34, No. 1, pp. 65-86. 

Maccagno, T. and Knott, J., 1991, “The Low Temperature Brittle Fracture Behavior of 
Steel in Mixed Modes I and II,” Engineering Fracture Mechanics, Vol. 38, No. 2/3, 

pp. 111 - 128 . 

Maccagno, T.M. and Knott, J.F., 1992, “The Mixed Mode I/II Fracture Behaviour of 
Lightly Tempered HY130 Steel at Room Temperature, ” Engineering Fracture Me- 
chanics, Vol. 41, No. 6, pp. 805-820. 

Maiti, S.K. and Mahanty, D.K., 1990, “Experimental and Finite Element Studies on 
Mode I and Mixed Mode (I and II) Stable Crack Growth-II. Finite Element Analy- 
sis,” Engineering Fracture Mechanics, Vol. 37, No. 6, pp. 1251-1275. 

Malvern, L.E., 1969, Introduction to the Mechanics of a Continuous Medium, Prentice 
Hall. 

Maurisch, T.D. and Ortiz, O., 1995, “A Finite Element Study of Chip Formation in High 
Speed Machining,” MED-Vol. 2-1/MH-Vol. 3-1, Manufacturing Science and Engi- 
neering, ASME, pp. 245-258. 

McMeeking, R.M. and Parks, D.M., 1979, “On Criteria for J-Dominance of Crack Tip 
Fields in Large-Scale Yielding,” ASTMSTP 668, American Society for Testing and 
Materials, Philadelphia, pp. 175-194, 1979. 

Nakagaki, M., Chen, W.H., and Atluri, S.N., 1979, “ A Finite-Element Analysis of Sta- 
ble Crack Growth - 1,” Elastic-Plastic Fracture, ASTMSTP 668, J.D. Landes, J.A. 
Begley and G.A. Clarke, Eds, American Society for Testing and Materials, pp. 195- 
213. ' 


88 



Newman, J.C., Jr., 1974, “Finite-Element Analysis of Fatigue Crack Propagation — In- 
cluding the Effects of Crack Closure,” Ph.D. Thesis, Virginia Polytechnic Institute 
and State University, Blacksburg, VA. 

Newman, J.C., Jr., Booth, B.C., and Shivakumar, K.N., 1988, “An Elastic-Plastic Finite- 
Element Analysis of the J-resistance Curve using a CTOD Criterion, ” Fracture Me- 
chanics: Eighteenth Symposium, ASTM STP 945, D.T. Read and R.P. Reed, Eds., 
American Society for Testing and Materials, Philadelphia, pp. 665-685. 

Newman, J.C., Jr., Dawicke, D.S., and Bigelow, C.A., 1992, “Finite-element Analyses 
and Fracture Simulation in Thin-Sheet Aluminum Alloy,” NASA TM 107662. 

Newman, J.C., Jr., Bigelow, C.A., and Shivakumar, K.N., 1993a, “Three-Dimensional 
Elastic-Plastic Finite-Element Analyses of Constraint Variations in Cracked Bodies, ” 
Engineering Fracture Mechanics, Vol. 46, No. 1, pp. 1-13. 

Newman, J.C., Jr., Dawicke, D.S., Sutton, M.A., and Bigelow, C.A., 1993b, “A Fracture 
Criterion for Widespread Cracking in Thin-Sheet Aluminum Alloys,” ICAF 17, In- 
ternational Committee on Aeronautical Fatigue, pp. 443-468. 

Ortiz, M. and Popov, E.P., 1985, “Accuracy and stability of integration algorithms for 
elastoplastic constitutive relations,” International Journal for Numerical Methods in 
Engineering, Vol. 21, No. 9, pp. 1561-1576. 

Ponthot, J.P. and Hogge, M., 1991, “The Use of the Eulerian-Lagrangian FEM in Metal 
Forming Applications Including Contact and Adaptive Mesh, "Advances in Finite 
Deformation Problems in Materials Processing and Structures, AMD- Vol. 1 25, 
ASME. 

Potyondy, D.O., 1993, “A Software Framework for Simulating Curvilinear Crack Growth 
in Pressurized Thin Shells.” Ph.D. Thesis, Cornell University, Ithaca, New York. 

Rice, J.R., 1968, “A Path Independent Integral and the Approximate Analysis of Strain 
Concentrations by Notches and Cracks,” ASME Journal of Applied Mechanics, Vol. 
35, pp. 379-386. 

Rice, J.R. and Rosengren, G.F., 1968, “Plane Strain Deformation near a Crack tip in a 
Power-Law Hardening Material,” Journal of the Mechanics and Physics of Solids, 
Vol. 16, pp. 13-31. 

Roy, S., Dexter, R.J., and Fossum, A.F., 1993, “A Computational Procedure for the 
Simulation of Ductile Fracture with Large Plastic Deformation,” Engineering Frac- 
ture Mechanics, Vol. 45, No. 2, pp. 277-293. 

Saouma, V.E. and Ingraffea, A.R., 1981, “Fracture Analysis of Discrete Cracking,” Col- 
loquium on Advanced Mechanics of Reinforced Concrete, International Association 
of Bridge and Structural engineers, Delft, June 1981, pp. 393-416. 

Schreyer, H.L., Kulak, R.F. and Kramer, J. M., 1979, “Accurate numerical solutions for 
elasto-plastic models f ASME Journal of Pressure Vessel Technology, 101, 226-234. 

Sekhon, G.S. and Chenot, J.L., 1993, “Numerical Simulation of Continuous Chip forma- 
tion During Non-Steady Orthogonal Cutting,” Engineering Computations, Vol. 10, 
pp. 31-48. 


89 


Shan, G.X., Kolednik, O., Stiiwe, H.P, and Fischer, D.F.,1992, “A Substitution Method 
for 3D Elastic-Plastic FE Analysis of Fracture Mechanics Specimens,” Engineering 
Fracture Mechanics, Vol. 41, No. 5, pp. 625-633. 

Shan, G.X., Kolednik, O., and Fischer, D.F., 1994, “A Numerical Study on the Influence 
of Geometry Variations on Stable Crack Growth in CT Specimens for Different Ma- 
terials,” Constraint Effects in Fracture: Theory and Applications, ASTMSTP 1244, 
Mark Kirk and Ad Bakker Eds., American Society for Testing and Materials, Phila- 
delphia. 

Shivakumar, K.N. and Newman, J.C., Jr., 1990, “ZIP3D - An Elastic-Plastic Finite Ele- 
ment Analysis Program for Cracked Bodies,” NASA TM-1 02753. 

Simo, J.C. and Taylor, R.L., 1985, “Consistent tangent operators for rate-independent 
elasto-plasticity,” Computer Methods in Applied Mechanics and Engineering, Vol. 
48, pp. 101-118. 

Sneddon, I.N., 1946, “The Distribution of Stress in the Neighborhood of a Crack in an 
Elastic Solid,” Proceedings, Royal Society of London, Vol. A- 187, pp. 229-260. 

Sommer, E. and Aurich, D., 1991, “On the Effect of Constraint on Ductile Fracture,” 
Proceedings of the European Symposium on Elastic-Plastic Fracture Mechanics, 
ESIS/EGF Pub. 9, pp. 141-175. 

Sutton, M.A., Zhao, W., Boone, M.L., Reynolds, A.P., and Dawicke, D.S., 1997, “Pre- 
diction of crack growth direction for mode I/II loading using small-scale yielding and 
void initiation/growth concepts,” International Journal of Fracture, No. 83, pp. 275- 
290. 

Swenson, D.V., 1985, “Modeling Mixed-Mode Dynamic Crack Propagation using Finite 
Elements,” Ph.D. Thesis, Cornell University, Ithaca, NY. 

Swift, T., 1993, “Widespread Fatigue Damage Monitoring Issues and Concerns,” 5th In- 
ternational Conference on Structural Airworthiness of New and Aging Aircraft.” 

Wawrzynek, P.A. and Ingraffea, A.R., 1987, “Interactive Finite element analysis of 
Fracture Processes: An integrated Approach,” Theoretical and Applied Fracture Me- 
chanics, No. 8, pp. 137-150. 

Wells, A. A., 1961, “Application of Fracture Mechanics at and Beyond General Yield- 
ing, ’’ British Welding Journal, Vol. 1 1, pp. 563-570. 

Westergaard, H.M., 1939, “Bearing Pressures and Cracks,” Journal of Applied Mechan- 
ics, Vol. 6. , pp. 49-53. 

Williams, M.L., 1957, “On the Stress Distribution at the base of a Stationary Crack, ’’ 
Journal of Applied Mechanics, Vol. 24, pp. 109-1 14. 

Yamada, T., 1993, “A Rezoning Procedure for Finite element Analysis of Rubber-like 
Materials,” Proceedings of the 5 th International Conference on Computing and 
Building Engineering, ASCE, pp. 327-334. 


90 




Equilibrium Equations 


In this section we consider the equilibrium of a general body and use the principle of 
virtual displacements to derive the finite element equations. The derivation follows from 
Malvern (1969) and Bathe (1996). Figure A.l shows the body under consideration. The 
displacement is prescribed on S u and tractions prescribed on S/. 



Figure A.l: A body in two-dimensions with prescribed displacements and tractions. 


The differential equations that describe the problem are equilibrium: 

°. (A - ,} 

throughout the body, with the natural (force) boundary conditions: 

<j n = f B (A.2) 

tj J Jl 9 

on S/, and the essential (displacement) boundary conditions: 
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We consider now any arbitrary, continuous virtual displacement Su t satisfying <5«, = 0 on 
S„. Then, 


+ = 0 (A.4) 

Since <5w, is arbitrary, (A.4) can be satisfied if and only if the term in parenthesis is zero. 
Using the chain rule: 


( CT jj$U, = O'gjSUj + CTySUjj , 

we have 

\y ~ + f! >SU ] dV = 0 - 

(A.5) 

Next, using the divergence theorem: 

l ),j = £ K<5w, ) n i dS , 

(A. 5) can be rewritten as 

l (-VySuv + fi H Su, )dv + £ (OySiij ) rij dS = 0. 

(A. 6) 

By introducing the boundary conditions (A.2) and (A.3), 

(<-<7 ,,Su,j + f‘Su,)dv + £ v /' Su-’dS = 0 . 

(A. 7) 

Also, because of symmetry of the stress tensor, (oj/=q/,), we have, 

= ^ J {H du i j +Su J =C7 >j S£ r 

(A.8) 

So the statement of virtual displacements becomes, 

l*„Se, i dV= | r f’& ll dV + ^f, s ’Su!’dS. 

(A.9) 



Equation (A.9) is valid for general equilibrium as long as a statically admissible stress 
distribution is defined, and as long as a kinematically admissible displacement solution is 
defined. The current work is limited to small strain, small rotation, but materially non- 
linear static equilibrium. In this case, the current volume corresponds to the initial vol- 
ume, and all stresses and strains are referred in the initial volume. The standard small 
strain measure is used: 

£ ,j = 2 ("/ j +»,,/)• 

We will now derive the governing finite element equations. We begin by discretizing 
(A.9) in the usual manner for finite elements (Bathe, 1996). First isoparametric interpo- 
lation functions for the displacements are defined as: 


u (m) =N (m) U 


(A. 10) 


where N (m) is the displacement interpolation matrix and U are the element nodal dis- 
placements for element m. The strains are interpolated using 


£<"'> = B (m) U 


(A.ll) 


where B (m) is the strain-displacement matrix obtained by the differentiating the dis- 
placement interpolation matrix, N (m) . 

We can write equation (A.9) in a discretized form for a single element, where the contri- 
bution of each element is added to the global system of equations: 


f Se {m)T <f (m) dV (m) = f Su sim ' T f s ' m 'dS lm) , 
Jk<"> is'"" 


(A. 12) 


where for brevity the body force term has been omitted. Substituting for the virtual 
strains and displacements, 


£u l [f ,. ) b ( " ) V" ) </f ( " ) }= 


(A. 13) 
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where SU is the variation of the displacement vector for the total domain. Since <5U is 
arbitrary (and dropping the superscript for the element, and still considering only one 
element): 

£B T <Ti/F= j^NfVS. (A. 14) 

A.1 Incremental Equations of Equilibrium 

The problems of interest here are commonly called quasi-static\ that is, there are no iner- 
tial effects. However, when nonlinear effects such as plasticity are considered, the solu- 
tion to the problem may be loading path dependent, and so, dependent on a pseudo-time. 

We introduce here the notion of time to represent equilibrium at discrete time points (eg. 
0, At, 2 At) where At is the increment in time. We assume that the full system state has 
been obtained up to time t, and we desire the solution at time t+At. Typically for quasi- 
static problems the time increment represents an increment in the loading conditions. 

A.1.1 Newton-Raphson Scheme 

One approach to solving non-linear finite element equations is the Newton-Raphson it- 
erative scheme. Following the development in Bathe (1996), we can rewrite (A. 14) as 

R(U*)= ,+A 'F e (U*)- ' +a 'F,(U*) = 0, (A. 15) 

where /+A 'F f (U*) is the external applied force vector, and ' +A, F I (U*) is the nodal internal 
force vector evaluated from the element stresses. At equilibrium R(U*) = 0. 

If we assume that the previous iterative solution ,+A, U <7-,) is available, then a Taylor series 
expansion gives: 

R(U’) = R(' +A 'U 0 -°) + — (U* - ,+A 'U (J - ,) ) + H.OT. (A. 1 6) 

9U 

Using (A. 15) we can rewrite (A. 16) as (neglecting higher-order terms): 

(U* _ _ ' +A 'F f _ ,+a 'F. u-I) (A. 17) 

The update equation is: 


3R 

SU 
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(A. 18) 


__/ + A/ j£(./-l)^^j(i) — /+/ ^ r p — /+ A/p(y-1) 

where /+A 'K (,_,) is the current tangent stiffness matrix: 

3F 




3U 




The iterative update to the displacement solution is 




(A. 19) 


(A.20) 


Several variations on the Newton iterative scheme exist. If (A. 19) is updated at each it- 
eration the scheme is typically referred to as the full Newton method; if (A. 1 9) is updated 
periodically, the scheme is typically referred to as the Modified Newton method. If only 
the initial value of (A. 19) is used, the scheme is typically referred to as the initial stress 
method. 

Efficiency concerns usually dictate which of the three is employed. The initial stress 
method requires many more iterations, but the Full-Newton and Modified-Newton meth- 
ods require that (A. 19) be evaluated and inverted more than once. A trade off usually 
exists between the convergence rate without updates to (A. 19), compared to the cost of 
reforming and inverting (A. 19). For small strain plasticity with no other nonlinearities, 
such as contact conditions, the initial stress method is usually adequate. The current im- 
plementation uses the initial stress method. 

A.2 Incremented Iterative Solution 

Accurate solutions of the nonlinear finite element equations rely not only on satisfaction 
of equilibrium, but also on accurate integration of the incremental elastic-plastic stress- 
strain constitutive relation. 

Crisfield (1991) provides a good description of implementing the stress update in terms 
of either incremental or iterative strains. Essentially, the choice lies in what reference 
state is used for the incremental update, and how the state variables are managed during 
the update. Crisfield references several approaches; only two are discussed here, with the 
most robust and conservative of the two being implemented. 
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Listing A. 1 shows the basic Newton steps for the iterative strain update method, where 
Sp is the iterative displacement computed from the inverse tangent stiffness matrix, 
K; 1 and the residual (out-of-balance) load vector g . The iterative strain update dz is 
calculated from the iterative displacement update, and the iterative stress update from the 
iterative strain update. All displacement strain, and stress variables reference the previ- 
ous iterative state (which is not in equilibrium). 

Listing A.2 shows the basic Newton steps for the incremental strain update method, 
where the iterative displacement dp is calculated as before. The algorithm departs here 
from the previous algorithm, in that the iterative displacement is accumulated into the 
incremental displacement, Ap, which is referred to the previous equilibrium state. 


1. 

8 . P = -K;'q 


2 . 

Sz = BAp 

= C,(<t) & 

3 . 

o„ = o 0 +da 


4 . 

Repeat 1-3 until convergence 


Listing A.1: Iterative strain update (Crisfield, 1991). 


Then the incremental strains, and incremental stresses can be calculated with respect to 

the previous equilibrium state. The new trial stress. 

<r„, is calculated from the old equilib- 

rium stress, ct 0 , and the stress increment referred to the previous equilibrium state, Ac . 
Once equilibrium is satisfied the trial stress is saved as the current stress. 

1. 

dp = -K”'q 

Ap = Ap + dp 

2 . 

Az - BAp 

A<t = C,(<r) Ae 

3 . 

<L, = «. + Ao 


4 . 

Repeat 1-3 until oynvergenoe 


5. 

<*o = 



Listing A.2: Incremental strain update (Crisfield, 1991). 


Figure 2 shows schematically for a one-dimensional problem a potential effect of using 
iterative vs. incremental updates. Point A represents a converged equilibrium state. Point 
B is determined from the first tangential iteration. In Figure 2(a) the second iteration pro- 
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duces point C, since the strain increment was negative, and Point B was used as a refer- 
ence. In contrast, in Figure 2(b), Point C is correctly determined by using the strain in- 
crement referred to Point A, so the negative strain iteration simply reduced the total strain 
increment. 

A.3 Implementation 

Bathe (1996) outlines nicely the steps necessary to implement the incremental Newton- 
Raphson solution scheme. Suppose that at time t equilibrium is established. The stresses 
'a, and strains ‘e, and internal material parameters 'k, are known. Using (A. 18), the initial 
displacement update (iteration /= 0), as 

£U (I) = ['K <,,) ]" l (' +A 'F e - ' +A, F/ n) ), (A. 21) 

where ,+A/ F e is the current external applied load and ' +a 'f/ 0) = ‘F^ n \ the converged internal 
forces from the previous load step at the final convergence iteration («). The displacent 
increments are updated using (A.20), and element strain increments are calculated using 

,+A 'A£ (i) = B ' +a 'AU (,) . (A.22) 

For elastic behavior the strains directly give the stresses using the appropriate stress- 
strain constitutive matrix. For inelastic behavior the stresses must be integrated as 

' +a V M) ='<i + | (A. 23) 

and the tangent stress-strain matrix ,+a 'D (M) must be evaluated consistent with the stress 
integration process. The stress integration, equation (A.21), is described in Appendix B. 


Once the stresses are updated, the nodal variables are updated using (A. 18) as 



(A.22) 

,+a 'AU' = ,+a, AU m + ,+A W . 

(A. 23) 

The total strains are evaluated as 


' +a, £ (/) = B ( ,+A 'U (0 + ,+a, AU (,) ). 

(A. 24) 
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The internal forces ' +A 'F, (i) are evaluated at the current trial stress state. Evaluating 


'+A/r(M) _ /+A/p(i) _ /-fA/p(/) (A. 25) 

e i 

ensures that any residual out-of-balance forces from the previous load step are included 
in the current iterations, thus preventing drift from the true equilibrium path, as might be 
experienced using a traditional initial stress formulation. 



Figure A.2: Schematic of (a) iterative and (b) incremental strain update (Crisfield, 1991). 
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Plasticity Rate Equations 


Appendix B 


The formulation of the plasticity rate equations follows in a similar manner to any of a 
number of finite element texts (Bathe, 1996; Crisfield, 1991). The plasticity rate equa- 
tions are formulated here for the von Mises yield criterion. The yield function for iso- 
tropic strain hardening is: 

f(< f > £ p*) = < T e- Voter,) ( Bl > 

where a is the current stress, s is the effective plastic strain cr e is the effective, or von 
Mises stress, and cr () is the yield stress. The strain rate is linearly decomposed into elas- 
tic and plastic parts: 

s =e e +i p ( B -2) 

where e e is the elastic strain rate and z p is the plastic strain rate. The flow rule is as- 
sumed to be associative: 

i=A^L = i a (B.3) 

p da 

The linear constitutive relation for the elastic stress rate is assumed to hold: 


ct = D(e — e p ) = D(e - Aa) 


(B.4) 


For isotropic strain hardening the equivalent plastic strain is 



(B.5) 


For plastic flow to occur (tangency condition) the stresses must remain on the yield sur- 
face: 


/ 


da dc 7 a de px ps 


(B6) 
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For the von Mises yield function under plane stress 


•2 


£ ps -J3 { £ px + £ py + £ px £ py 4 Y pxy J 


r „ - (o ', 2 + <V - <r, <T, + IcrJ Y* 

if f 

4 / pxy ) 

Now the flow rule becomes 


K -k s ±- 

p da 


(B.7) 

(B.8) 



i 

2(7 „ - <J r - 

€ n%) 


py 

2 G e 

y x 


e 

6cr 

1 m J 


*y J 


(B.9) 


Under uniaxial tension, a x , s ny = £ p2 =~'A£ rx , so that there is no plastic volume change; 

= <j x . Cons 

be taken from the uniaxial stress / plastic strain relationship. 


so, £ p =£ px , and cr e = a 0 = a x . Consequently, the relationship between cr 0 and s ps can 


H > = <*Ll. _ - JL 


te px ten* 1 " 7e 


(B.10) 


ps 


Now the tangency condition is 


• df\ 8f J dcr 0 . t • tji- n 
J da dcr n d£„ p ' p 


(B.l 1) 


ps 


Substituting (B.9) into (B.8) gives 


£ P s = ^ 


(B.12) 


and substituting (B.12) into (B.l 1) gives 

/ = a 7 a - A'X 


(B.13) 


where A' = A\s ) is substituted for H' and is in general a non-linear hardening 
modulus that is a function of the plastic strain. Premultiplying (B.4) by a T and substi- 
tuting into (B.13) gives 


X = 


a T De 


a T Da + A' 

and taking into account symmetry in D, 


(B.14) 
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(B.15) 


f 


* = D„e = D| 


aa T P 
a T Da + A' 


\ 

8 

J 


This formulation cannot account for the Bauschinger effect induced by cyclic loading. 
The initial analyses under consideration require only monotonic loading, so the Bausch- 
inger effect will not be introduced. A kinematic hardening model may be introduced at 
some later time to account for the ‘kinematic shift’ of the yield surface in this case. 


B.1 Integrating the Plasticity Rate Equations 

The analyses considered here are for quasi-static conditions and consider no dynamic ef- 
fects, so the dotted quantities can be considered small changes. The time rate of change 
relations can then be written in the form de = edt . The other dotted quantities can be 
considered similarly. For this section, the rate equations will be considered in this incre- 
mental form, rather than rate form. 


Many algorithms are available for integrating the incremental form of the von Mises con- 
stitutive equations. Many are based on what has become known as the generalized trape- 
zoidal rule (Ortiz and Popov, 1985; Crisfield, 1991). From the generalized trapezoidal 
rule can be derived the pure ‘explicit’ or ‘forward-Euler’ scheme, a mid-point algorithm, 
or an ‘implicit’ or ‘backward-Euler’ scheme. Accuracy of the integration algorithm is of 
prime concern, with efficiency also a concern. Krieg and Krieg (1977), Schreyer, et. al. 
(1979), Ortiz and Popov (1985), and De Borst and Feenstra (1990) all showed that for 
large step sizes the backward-Euler scheme is most accurate, and in certain cases is by far 
the most efficient, since for plane strain with linear hardening the backwards-Euler 
scheme degenerates into the well know radial return algorithm (Wilkins, 1964; Krieg and 
Key, 1976). The basic algorithm implemented uses an elastic predictor, normal corrector. 
Figure B.l is a schematic of the elastic predictor, normal corrector implemented. The 
scheme uses the normal at the ‘elastic trial point,’ B. A first order expansion about point 
B of the yield function gives: 

/ = / B + 4* + AL A £ = /„ - AA*lC*. - AAA’ (B. 1 6) 

So d£ ps 

Solving for A/l 


A/l = 


Jb 

a T B Da B + A' h 


Then the stresses at C are given by 


(B. 17) 
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(B.18) 


o c = a B - AADa B 

For the 3-D (or plane strain or axisymmetric) von Mises with linear hardening, (B.18) 
corresponds to the radial-return algorithm, which is a special form of the backward-Euler 
algorithm. For non-linear hardening, or for plane stress, the stress at o c will not neces- 
sarily lie on the yield surface, and iterations are required to satisfy the consistency condi- 
tion. 

The backward-Euler return implemented here comes from Crisfield (1991). A simple 
statement of the equation, referring to Figure B.l is: 

<r c = <r B - A/lDa c (B.l 9) 

where o B is the elastic trial stress and the initial value for a c comes from the value in 
(B.18). An iterative method is required to solve equation (B.l 9), since the normal at the 
trial position B will not necessarily equal the final normal. A Newton type iterative 
scheme is implemented, where the residual vector of stresses, r, is written as: 

r = < 7 -(a B - ATDa c ) (B.20) 

and iterations reduce the residual to (almost) zero. The final stresses a should satisfy the 
yield criterion, /= 0. 
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Figure B.1: Elastic predictor, normal return with correction to the yield surface for plane stress. 


With the trial stresses ct b kept fixed, a truncated Taylor series expansion is applied to 
produce a new residual, r n as: 

r = r + d + ADa + AAD — d (B.2 1 ) 

n 0 do 


where d is the change in a and X is the change in AA . Setting r n to zero gives 


o = - 


I + AID (r 0 + ADa) = -Q r 0 - iQ“'Da 


(B.22) 


Also, a truncated series expansion on the yield function gives 

f = f + — d +-^—e n . = / + aid - A'X = 0 
Jn Jco do de ps ps c ° c 

Solving for i from (A 22) and (A 23) and dropping the subscript C gives 
; /„-a T Q-'r„ 

f — /l — T , 

PS ^ T -1 1 


(B.23) 


a Q Da + A' 


(B.24) 


Substituting (B.24) into (B.22) gives the iterative change in the stress, and (B.24) is the 
iterative change in effective plastic strain. It is worth noting here that at each iteration the 
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yield function is evaluated at /„ = f„( s PH + & £ P .s) and the hardening at 
A' = A'(e pSB + Ae p J, where As p!i is the accumulated change in e ps from point B. The 
term dsi! ct t is conveniently derived as follows from the yield equation: 



(B.25) 

v =a= i 

(B.26) 

da 2a e da 


<*> _ i i aa T 

(B.27) 

da 2cr e da 2 a e 



B.2 Consistent Tangent Modulus 

The current implementation of the equilibrium iterations uses the initial stiffness matrix 
in the modified Newton-Raphson scheme. If a full Newton scheme is required for accel- 
erated convergence a tangent modular matrix that is fully consistent with the backward- 
Euler integration algorithm can significantly improve convergence. Simo and Taylor 
(1985) derived a tangent modular matrix that is fully consistent with the backward-Euler 
integration algorithm. Standard techniques would use the modular matrix of (B.15), 
which is inconsistent with the backward-Euler algorithm, and as such, does not provide 
the quadratic convergence associated with the Newton-Raphson method. More recently 
Crisfield (1991) derived a consistent tangent modular matrix based on the above equa- 
tions for integrating the plasticity rate equations. This derivation begins with the expres- 
sion for the backward-Euler algorithm: 

<t = o B - ATDa (B.28) 

where we drop the suffix C for the current configuration, but the suffix B is kept to indi- 
cate the elastic trial stress. Differentiation of (B.28) gives 

<t = Dd - TDa - AAD — a (B.29) 

da 

The last term of equation (B.29) is omitted from the derivation of the standard tangent 
modular matrix in (B. 1 5). Now, from (B.29) we have 

a = (l + AAD— ) d(e - Aa) = Q"'d(s - ia) = r(e - ia) (B.30) 
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where the Q matrix appeared before in (B.22). 

For f to remain on the yield surface it’s value should be zero, so from (B.6) and (B.30) 
f = a T d - A'A = a T Re - ia T Ra = 0 (B.3 1 ) 


Solving for k and substituting into (B.30) gives 


d = D ct e = 


T r> T 


R- 


Raa R 


a T Ra + A') 

where account of the symmetry of R has been taken in deriving equations (B.32). 


(B.32) 
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Appendix C 


Plasticity Verification 


This example demonstrates a number of FRANC2D/L features for a materially-nonlinear 
analysis, including the von Mises material model with linear isotropic hardening and 
nonproportional loading. A thin tube is first stressed in tension to the point of yielding 
and is then twisted under constant axial stress. A single material point is sufficient to 
model the problem. Figure C.l shows the loading path. 


x 


40 


O 



Point 

a 

T 

O 

0.0 

0.0 

A 

26.25 

0.0 

B 

26.25 

24.0 


Figure C.l: Loading path for the tension torsion example. 


For linear isotropic hardening the Hill (1950) plasticity book gives the exact solution for 
the axial and shear strains as: 


V 

f 3r 2> l 

Y„ 

0 

-In 1 + — H 

+ -f- 

2 H' v y; j 

E 

3 

Y a 



t--i= tan 
V3 

l n Jj 


where Y 0 is the yield stress in tension, E is Young 
and G is the elastic shear modulus for v= 0.3. 


(C.l) 

(C.2) 

’s Modulus, H’ is the plastic modulus, 


Since a single material point is sufficient to model the problem, only one Q8 element in a 
single layer is necessary for a FRANC2D/L model. Figure C.2 shows the material prop- 
erties used, and Figure C.3 shows the applied boundary conditions. The loads were ap- 
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plied in two separate load cases. Then load factors were used to scale the loads. The ten- 
sile load factor was set to 26.25 and the shear load factor to 0.0, then equilibrium was es- 
tablished. Next the tensile load factor was set to 0.0 (load held constant, no increment) 
and the shear load factor was slowly incremented to apply the twist under constant ten- 
sion. 



Figure C.2: Material properties for the tension torsion example. 


Figure C.4 shows the exact strain path and FRANC2D/L results for selected applied 
loads. The results show that the FRANC2D/L plasticity model is capable of accurate re- 
sults. 


A O O 


CT 



T 


•+ 


a 


a = 1.0 
x — 1.0 


Figure C.3: Fixity and loading conditions. Shear loads and axial loads are in different load cases 
to facilitate non-proportional loading. 
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Elastic-Plastic Validation 


Appendix D 


Modeling elastic-plastic fracture requires an accurate description of strains in the process 
zone. Newman et al. (1992) compared elastic-plastic finite element analysis results with 
experimental results for a blunt-notch specimen. These specimens were similar to M(T) 
specimens, except that a small hole was drilled at the end of the saw cut, and the holes 
polished to help prevent premature fracture. 

Figure D. 1 shows the geometry of the specimen, and Figure D.2 shows the finite element 
idealization. Only !4 of the specimen was modeled by using symmetry along the hori- 
zontal and vertical center-lines of the specimen. Displacement control boundary condi- 
tions were applied along the top edge of the model. Table D.l contains the analysis ma- 
terial properties for 2024-T3 aluminum, and Table D.2 contains the multi-linear stress- 
strain curve used in the von Mises yield model. 


Table D.l: Analysis Material Properties 


Young’s Modulus 

10350 Ksi 

Poisson’s Ratio 

0.3 

Thickness 

0.09 inch 

Yield Stress 

50 Ksi 


Figure D.3 is a contour plot of the von Mises stress near the failure load. The plot shows 
the extent of yielding in the model. Figure D.4 shows the applied stress against notch- 
root displacement for the test and for the FRANC2D/L simulation. Displacements were 
measured at both notch roots using displacement gages. The FRANC2D/L results show 
excellent agreement with the experimental measurements, indicating that the small strain 
software is capable of simulating the relatively large-scale plasticity experienced by the 
test specimen. 
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Figure D.1: Schematic of blunt notch specimen (w = 10 in, h = 12 in, c = 1.875 in). 



Figure D.2: Finite element idealization of blunt notch specimen. 
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Table D.2: Uniaxial stress-strain curve. 


Stress (Ksi) 

Strain 

0.0 

0.0 

50.0 

0.00483 

56.5 

0.015 

62.5 

0.04 

68.5 

0.1 

71.0 

0.16 

1.0 

0.2 


1 71.30 
64.02 
56 73 
49.45 
42 17 
34 68 
27.60 
20 32 
13 03 
5.751 
- 1.532 




Figure D.3: Contour plot of von Mises stress (ksi). 
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Figure D.4: Blunt notch configuration and results. 
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